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ABSTRACT 


This  study  analyzes  the  statistics  of  sensor  data  obtained  from  the  ob¬ 
servations  of  artificial  satellites.  The  effect  of  these  statistics  on  the 
estimation  of  the  orbital  elements  has  been  evaluated. 

Several  techniques  for  removal  of  the  actual  orbit  from  the  observa¬ 
tions  to  derive  residuals  have  been  investigated.  The  residuals  have 
been  obtained  using  both  polynomial  fitting  with  and  without  autoregressive 
scheme.  From  the  satellites  considered,  no  consistency  was  apparent 
in  the  magnitude  of  the  variances  of  the  individual  parameters.  How¬ 
ever,  a  consistency  was  noted  in  the  skewness  of  the  residual.  Com¬ 
puter  programs  have  been  developed  for  the  determination  of  the  statis¬ 
tics  of  the  orbital  elements  from  observation  data.  The  study  shows 
that  the  variances  of  the  orbital  elements  are  functions  of  the  variance 
of  the  observations,  the  time  of  observation  and  description  of  the  mo¬ 
tion  of  the  satellite.  Furthermore  the  variance  of  any  orbital 
elements  is  minimum  if  the  sum  of  the  variance  of  the  observations  at 
a  given  time  is  minimum  for  a  specific  computational  procedure  of  the 
orbital  element  from  observations. 

The  behavior  in  the  statistics  of  the  orbital  elements  exhibits  a  consist¬ 
ency  in  the  measure  of  skewness  and  excess.  Since  these  quantities  are 
normalized  with  respect  to  the  standard  deviations  it  appears  that  a  bias 
is  introduced  in  the  computation  procedure.  Thus,  the  problem  of  ob¬ 
taining  accurate  position  prediction  of  artificial  satellites  is  one  of 
interplafy  between  the  quality  of  observations  and  an  analytical  description 
of  the  motion  of  the  satellite.  How  to  trade  these  factors  depends  upon 
the  time  allowed  for  deriving  the  orbit  from  observations. 
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INTRODUCTION 


This  final  report  summarizes  the  results  of  a  study  on  the  error  of 
orbital  estimation  from  observed  data  performed  by  the  Radio  Corpora¬ 
tion  of  America  for  SPO  496 L  Space  Track  R&O  facility  under  Contract 
AF19(604)  -8041. 

Determining  accurate  orbits  for  an  artificial  satellite  has  presented  a 
new  problem  in  celestial  mechanics.  To  imply  that  accurate  positional 
information  of  an  orbiting  body  can  be  obtained  by  more  accurate  obser¬ 
vation  is  only  partially  correct.  There  is  need  for  an  analytical  theory 
that  would  represent  the  motion  of  the  vehicle  %vith  an  accuracy  consis¬ 
tent  with  the  accuracy  of  the  observations. 

This  study  was  concerned  mainly  with  the  accuracy  of  the  observations 
and  their  effects  on  the  orbital  elements.  More  specifically,  the  study 
was  devoted  to 

(1)  The  analysis  of  the  statistics  of  satellite  tracking  observations 
using  data  from  tracking  and  surveillance  radars  and  Baker- 
Nunn  cameras, 

(  2)  The  investigations  of  the  statistical  and  dynamical  assumptions 
governing  the  orbit  computation  based  on  the  Herrick-Gibbs 
method  and  differential  correction, 

(3)  The  analysis  of  the  error  propagation  from  observations  to 
orbital  elements  following  the  Herrick-Gibbs  method, 

(4)  The  investigation  of  intermixing  different  types  of  data  and 
criteria. 
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Determination  of  the  statistical  quality  of  satellite  tracking  observations 
was  conducted  by  gathering  satellite  tracking  data  essentially  from  the 
Millstone  Hill  Radar  and  AN/FPS-49.  This  analysis  is  contingent  on 
the  removal  of  the  actual  orbit  from  the  observations.  Several  tech¬ 
niques  were  employed.  Section  1  of  this  report  discusses  these  methods; 
namely,  simulation,  differential  correction,  and  polynomial  smoothing 
with  and  without  an  autoregressive  scheme.  Each  of  these  methods  re¬ 
quires  great  care  in  its  application,  and  the  statistics  gathered  should 
not  be  regarded  as  absolute. 

The  orbit  removal  technique  depends  upon  assumptions.  To  determine 
whether  the  assumptions  adequately  describe  nature,  it  is  necessary  to 
look  at  the  results  of  the  orbit  removal- statistical  analysis  procedure. 
There  are  two  types  of  tests  that  may  be  applied.  First,  it  may  be 
noted  whether  the  statistics  are  self- consistent, i.  e. ,  is  a  different 
statistical  behavior  for  a  given  sensor's  observations  obtained  on  dif¬ 
ferent  satellites? 

The  second  test  is  the  "how  well  does  it  work?"  type.  Since  the  end 
product  of  the  statistical  analysis  is  to  be  a  scheme  for  precision  deter¬ 
mination,  estimates  of  sensor  statistics  may  be  used  to  construct  an 
"optimum  estimator,  "  and  to  test  the  orbit  so  estimated  against  what¬ 
ever  high  quality  observations  are  available. 

A  large  number  of  restrictions  are  placed  on  the  use  of  differential 
correction.  Without  drag  correction  the  observational  span  must  be 
restricted  to  a  few  days,  and,  in  any  case,  to  a  few  weeks.  A  very 
limited  number  of  observations  from  each  sensor  must  be  used  and 
the  relative  quality  of  observations  from  different  sensors  must  be 
known.  Furthermore,  a  fairly  large  number  of  observations  is  neces¬ 
sary  in  order  to  obtain  statistically  stable  residuals. 
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These  restrictions,  of  course,  greatly  limit  the  number  of  satellites 
to  which  this  technique  can  be  applied.  However,  there  will  probably 
be  some  satellites  with  sufficient  concentration  of  data.  Fox  these 
satellites,  the  amount  of  computation  necessary  is  large,  especially  if 
the  iterative  procedure,  suggested  in  subsection  1.2.  1,  is  applied.  An 
alternative  approach  to  the  orbit  removal  problem  is  polynomial  smooth¬ 
ing  which  requires  much  less  computation. 

The  simpler  method  of  polynomial  smoothing  requires  concentrated 
data.  This  method,  which  has  been  applied  to  Millstone  data,  has  been 
reasonably  successful.  However,  two  procedures  must  be  noted.  In 
the  first  place,  the  observations  to  be  reduced  have  to  be  taken  over  a 
time  interval  limited  by  the  accuracy  of  the  satellite  theory.  For  ex¬ 
ample,  if  it  is  desired  to  analyze  range  observations  from  Millstone 
which  have  a  standard  deviation  on  the  order  of  one  kilometer,  then 
these  observations,  plus  any  others  which  are  used  to  establish  the 
orbit,  must  cover  a  time  span  less  than  that  over  which  the  satellite 
theory  is  accurate  to,  say,  one-half  kilometer.  This  time  span,  even 
for  a  high  satellite,  is  probably  not  much  more  than  a  few  days  at  the 
most. 

Secondly,  using  observations  from  different  sensors  without  properly 
weighting  them  should  be  avoided.  Since  the  statistical  analysis  to  be 
performed  is  a  prerequisite  to  the  determination  of  weights,  it  seems 
best  to  perform  the  correction  on  observations  made  by  a  single  sensor. 

A  disadvantage  of  the  polynomial  fitting  is  that,  in  the  topocentric  carte¬ 
sian  coordinates,  all  the  parameters  are  smoothed  by  a  same  degree 
polynomial.  The  transformation  to  cartesian  coordinates  is  an  attempt 
to  obtain  uncorrelated  residuals.  With  the  autoregressive  scheme  em- 
polyed  in  conjunction  with  the  polynomial,  the  smoothing  is  done  directly 
on  the  radar  parameters.  Each  one  can  then  be  fitted  by  a  different 
optimum  degree  of  polynomial.  The  satellite  observations  handled  by 
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this  method  indicate  that  the  residuals  of  the  radar  parameters  were 
uncorrelated. 

In  the  process  of  this  investigation  it  was  found  that  the  differential  cor¬ 
rection  procedure  used  for  astrodynamical  data  reduction  is  to  be  re¬ 
garded  as  a  method  of  solution  of  the  unweighted  least  square  problem 
(Ref.  2).  It  was  also  found  that  the  statistics  of  the  residuals  exhibit 
non-zero  third  and  fourth  moments.  It  can  be  shown,  using  the  fitting 
of  observations  by  a  trend  and  autoregressive  scheme,  that  best  esti¬ 
mators  for  predicting  orbital  elements  can  be  obtained  with  independent 
estimators  of  the  observations  (Ref.  10).  Expressions  to  obtain  the 
statistics  of  the  osculating  elements  and  error  propagation  have  been 
derived  analytically.  Their  usefulness  is  limited  due  to  the  cumbrous 
formula.  A  computer  program  has  been  devised  to  derive  the  statistics 
of  the  osculating  elements.  This  program  utilizes  the  statistics  of  the 
observations,  and  the  Herrick-Gibbs  method.  It  gives  the  moments  of 
the  osculating  elements. 
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SECTION  1 


SMOOTHING  AND  DATA  ANALYSIS 

To  determine  the  statistical  character  of  electronic  or  optical  observa¬ 
tions  of  artificial  earth  satellites,  it  is  necessary  to  isolate  the  observa¬ 
tional  errors  from  the  "ideal”. observations.  The  observations  may  be 
of  slant  range,  azimuth  and  elevation,  or  of  right  ascension  and  declina¬ 
tion.  Therefore,  one  must  determine,  implicitly  or  explicitly,  the  true 
orbit  of  the  satellite,  calculate  from  this  orbit  the  values  which  each 
sensor  should  have  observed  (the  ideal  observations),  and  subtract  the 
ideal  from  the  actual  observations  to  obtain  the  errors. 

The  ideal  observations  may  be  determined  by  standard  astrodynamical 
methods  if  the  satellite's  orbit  is  known.  However,  in  order  to  use  an 
orbit  to  obtain  accurate  observational  errors,  it  is  necessary  to  know 
the  orbit  with  extremely  great  precision.  Since  the  orbit  is  usually  de¬ 
termined  by  applying  a  differential  correction  reduction  on  various 
observations,  certain  procedures  must  be  avoided  in  this  reduction  or 
the  "observational  errors"  obtained  become  meaningless. 

1.  1  DETERMINATION  OF  PARAMETRIC  STATISTICS 

1.1.1  DESCRIPTION  OF  METHODS 

The  inherent  statistics  associated  with  any  observation  parameter  ob¬ 
tained  from  the  sensing  equipment  should  be  fully  understood,  together 
with  a  knowledge  of  their  magnitudes,  limitations  and  model  assump¬ 
tions  presumed. 
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These  statistics  could  be  obtained  either  from  a  data  analysis  study  of 
a  simulated  target  or  from  actual  satellite  observations.  Both  techniques, 
unfortunately,  have  their  own  advantages  and  disadvantages. 

1.1.2  SIMULATION  STUDY 

Parameter  statistics  are  sometimes  gathered  by  means  of  simulated 
targets.  For  electronic  sensors,  signals  representing  those  expected 
from  an  actual  satellite  may  be  introduced  at  the  front  end  of  the  re¬ 
ceiver.  When  this  is  done,  the  orbit  removal  problem  becomes  trivial, 
for  the  ideal  observations  are  determined  by  the  "orbit"  which  was 
simulated. 

Besides  reducing  the  problem  of  orbital  removal  to  a  purely  numerical 
one,  the  technique  of  target  simulation  has  several  other  important  ad¬ 
vantages.  As  much  statistical  data  as  desired  may  be  gathered,  merely 
by  replaying  the  target  through  the  sensor- computer  system.  Further, 
the  quality  of  the  observations  may  be  calculated  as  a  function  of  signal 
strength.  Another  important  characteristic  is  that  simulation  offers 
a  powerful  tool  for  separating  the  random  and  systematic  components 
of  error. 

The  chief  disadvantage  of  simulation  is  that  one  cannot  be  sure  all 
sources  of  error  are  included.  Propagation  errors,  ground  reflections 
and  target  scintillation  are  usually  ignored.  Furthermore,  special 
equipment  used  for  generating  the  signal  may  introduce  unwanted  errors. 
In  spite  of  these  limitations,  it  is  felt  that  statistical  results  gathered 
from  a  sophisticated  simulation  program  should  provide  enough  infor¬ 
mation  on  the  character  of  observational  errors  to  greatly  simplify  the 
analysis  of  tracking  data.  Unfortunately,  extensive  simulation  has 
been  carried  out  for  very  few  sensors. 
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1.1.3  TRACKING  DATA  ANALYSIS 


To  obtain  meaningful  statistical  information  from  actual  satellite  obser¬ 
vations,  a  large  number  of  observations  must  be  available.  However, 
because  of  the  errors,  no  possible  satellite  orbit  will  exactly  satisfy  all 
of  the  observations.  Thus,  in  order  to  estimate  the  elements  of  the 
orbit,  some  criterion  is  needed  for  closeness  of  fit  of  actual  and  ideal 
observations.  This  criterion  is  related  to  the  mathematical  model 
adopted  for  the  errors,  in  the  sense  that  different  models  give  rise  to 
different  criteria  for  efficient  estimators  of  the  elements,  and  the  dif¬ 
ferent  estimators  will  usually  give  rise  to  different  sets  of  observational 
residuals.  If  an  error  model  which  does  not  adequately  describe  the 
character  of  the  errors  is  used,  the  residuals  obtained  will  be  signifi¬ 
cantly  at  variance  with  both  the  adopted  and  the  true  model.  The  analy¬ 
sis  of  satellite  tracking  data  may  be  regarded  as  an  iterative  procedure. 
In  this  procedure  an  orbit  is  determined  which  minimizes  a  criterion 
based  on  the  assumed  error  model  and  the  residuals  are  tested  against 
the  model.  A  good  first  approximation  is  needed,  since  "convergence" 
is  rarely  assured. 

1.2  DIFFERENTIAL  CORRECTION 

1 . 2.  1  DESCRIPTION  OF  MODEL  FOR  DIFFERENTIAL  CORRECTION 

The  correct  selection  of  an  error  model  to  best  describe  the  character 
of  the  errors  is  extremely  important  if  the  residuals  are  to  have  signifi¬ 
cant  meaning.  The  simplest  error  model  consists  of  the  following 
assumptions: 

(1)  All  observational  errors  are  Gaussian  random  variables 

(2)  All  errors  are  uncorrelated 

(3)  All  errors  have  zero  means  (observations  are  unbiased). 
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With  this  model,  the  most  "efficient''  estimator  is  the  rigorous  least 
squares  estimator.  An  exacting  discussion  on  the  most  "efficient" 
estimator  is  given  in  Ref.  2.  The  estimator  picks  out  that  orbit  which 
minimizes  the  sum  of  the  squares  of  the  weighted  residuals,  where  the 
weight  for  a  given  residual  is  the  inverse  of  the  standard  deviation  of 
the  observational  error.  To  actually  compute  the  orbital  element  esti¬ 
mates,  a  nonlinear  system  of  equations  must  be  solved.  Differential 
correction  is  an  iterative  technique  for  solving  these  equations. (See 
Ref.  1 ,  2  and  3. ) 

As  mentioned, the  weights  to  be  used  with  a  differential  correction  de¬ 
pend  iqxxi  the  standard  deviations  of  the  errors.  Since  these  standard 
deviations  are  to  be  obtained  by  the  statistical  analysis  which  follows 
the  orbit  removal,  the  need  of  an  additional  assumption  describing  the 
relative  sizes  of  the  error  is  apparent.  This  assumption  is  crucial, 
and  will  be  discussed  in  detail  in  subsection  1. 2.  3. 

1. 2.  2  THEORETICAL  LIMITATIONS  OF  DIFFERENTIAL  CORRECTION 

To  obtain  meaningful  residuals,  the  differential  correction  procedure 
must  not  introduce  theoretical  errors,  i,  e. ,  the  correction  must  select 
the  true  least  squares  estimate  of  the  orbit  to  rather  high  accuracy.  In 
particular,  this  requirement  imposes  the  following  restrictions  on  the 
procudure: 

(1)  Quantities  derived  from  observations  must  be  corrected  if 
their  uncertainty  gives  rise  to  significant  variance  in  the 
residuals 

(2)  The  time  span  of  the  observations  must  be  no  greater  than  the 
span  over  which  the  satellite  theory  is  of  "sufficient  accuracy" 

(  3)  The  loss  of  accuracy  involved  in  computing  partial  derivatives 
from  a  Keplerian  ellipse  assumption  must  be  investigated 
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(4)  Great  care  should  be  taken  in  the  numerical  processes,  especial 
ly  the  matrix  inversion. 

The  first  restriction  can  best  be  explained  by  a  numerical  example.  One 
effect  of  atmospheric  drag  is  to  reduce  the  period  of  a  satellite.  A  typi- 
cal  value  of  this  reduction  is  10  days  per  revolution  per  revolution. 
Suppose  a  preliminary  estimate  of  the  drag  parameter  has  been  made 
which  is  sufficiently  uncertain  as  to  make  the  value  of  the  period  reduc¬ 
tion  uncertain  up  to  10“^  days/rev^.  This,  in  turn,  causes  a  variable 
time  uncertainty  whose  magnitude  is  given  approximately  by  10  times 
the  square  of  the  number  of  revolutions  from  the  center  time.  For  two 
weeks  of  observations  this  amounts  to  a  maximum  time  error  of  about 

_5 

10  days,  hence  a  maximum  position  error  of  about  6  or  7  kilometers. 

Such  an  error  is  larger  than  the  standard  deviation  of  most  observations; 

therefore,  it  is  necessary  to  correct  the  drag  parameter  as  well  as  the 

-9 

orbital  elements.  Indeed,  with  an  uncertainty  of  10  drag  corrections 
are  necessary  whenever  the  observations  span  more  than  about  three 
days,  or  even  less  time  if  high-accuracy  optical  observations  are  to  be 
analyzed. 

The  second  restriction  refers  to  the  state  of  the  art.  If  exact  initial 
conditions  of  the  satellite  at  some  epoch  are  given  and  the  position  of 
the  satellite  can  be  computed  to  an  accuracy  of  0.  Z  km  for  only  one 
week  in  the  future,  then  the  observational  span  should  be  restricted  to 
about  two  weeks  or  less.  This  restriction,  unlike  the  other  three, 
cannot  be  removed  by  improving  the  differential  correction  program. 

The  third  restriction  is,  perhaps,  less  important  than  the  others.  The 
differential  correction  procedure  employs  the  partial  derivatives  of 
ideal  observations  with  respect  to  the  parameters  which  are  to  be  cor¬ 
rected.  These  derivatives  may  be  computed  either  numerically  or 
analytically.  When  analytical  differentiation  is  employed,  the 
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derivatives  are  computed  on  the  basis  of  a  two-body  Keplerian  orbit. 

This  introduces  a  second  order  error  into  the  estimates.  The  fourth 
restriction  is  self-explanatory. 

1.2.3  STATISTICAL  LIMITATIONS  OF  DIFFERENTIAL  CORRECTION 

The  properties  which  determine  the  degree  of  sophistication  necessary 
in  the  differential  correction  program  have  been  examined.  The  type 
of  data  to  be  used  with  the  program  is  analyzed  in  the  light  of  the  as¬ 
sumptions  on  the  error  model. 

The  assumption  of  Gaussian  error  distribution  is  not  too  critical.  Actual 
distributions  of  errors  will  not  differ  too  much  from  the  Gaussian,  and 
for  these  distributions  the  least  squares  estimator  will  still  be  highly 
efficient.  The  assumptions  of  unbiased  and  statistically  uncorrelated 
observations  are  more  important.  If  these  assumptions  are  not  satis¬ 
fied,  and  in  practice  they  may  not  be,  the  accuracy  of  the  estimated 
orbit  is  greatly  reduced.  For  this  reason,  whenever  there  is  a  sus¬ 
picion  of  such  bias  or  correlation,  the  amount  of  data  used  from  each 
pass  over  each  sensor  should  be  restricted  to  a  few  points.  In  this 
way,  the  bias  may  be  treated  as  an  addition  to  the  random  error,  and 
the  effects  of  serial  correlation  are  greatly  reduced. 

The  assumption  on  the  relative  quality  of  observations  is  critical. 

Briefly,  the  accuracy  of  an  improperly  weighted  least  squares  estima¬ 
tion  is  determined  by  the  more  grossly  overrated  observations.  Thus, 
if  the  data  include  several  observations  of  very  poor  quality,  these 
observations  must  have  corresponding  low  weights  attached  or  they 
seriously  distort  the  residuals. 

For  this  reason,  if  good  a  priori  information  on  the  relative  quality  of 
the  observations  is  not  available,  the  differential  correction  should  be 
regarded  as  an  iterative  procedure,  whereby  the  residuals  are  examined 
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to  determine  new  weights  for  the  next  correction.  Such  a  procedure 
has  no  guarantee  of  convergence. 

1.3  POLYNOMIAL  SMOOTHING 

The  technique  of  polynomial  smoothing  can  be  used  when  compact  sets 
of  satellite  observations  are  available.  The  polynomial  smoothing  con¬ 
sists  of  finding  those  time  polynomials  of  a  specified  degree  which  best 
fit  the  azimuth,  elevation  and  range  parameters  in  the  least  square 
sense.  The  polynomial  smoothing  technique  is  based  on  the  same  as¬ 
sumptions  that  apply  to  differential  correction.  The  former  technique 
gives  a  less  optimum  estimate  but  requires  much  less  computation 
than  the  latter.  Polynomial  smoothing  is  applicable  only  to  radar 
tracking  data,  whereas  the  differential  correction  technique  may  be 
applied  to  a  wide  variety  of  observations.  The  observations  need  not 
be  taken  at  equal  intervals  of  time. 

1. 3.  1  DESCRIPTION  OF  POLYNOMIAL  SMOOTHING 

The  polynomial  smoothing  technique  consists  of  finding  that  poly¬ 
nomial  of  a  given  degree  which  best  fits  radar  observations  of  a  given 
type  in  the  least  squares  sense.  For  example,  to  fit  radar  range  with 
a  second  degree  polynomial,  three  quantities  _  Cq,  Cj,  _  must 

be  found  which  minimize  the  expression,  S: 

S(C(„Cj.C3).^(R^(t)-Co-Cj  t.c/)2 

t 

where  is  the  observed  range  at  time,  t;  the  inverse  of 

the  standard  deviation;  and  the  summation  extends  over  all  the  obser¬ 
vation  times.  This  problem  has  a  unique  solution  unless,  of  course, 
there  are  fewer  observations  than  constants  to  be  determined.  The 
solution  reduces  to  a  linear  system  of  equations.  The  application  of 
this  technique  to  the  problem  of  orbit  removal  is  presently  discussed. 
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The  range  observations  are  the  sums  of  the  ideal  observations  plus 
the  observational  errors.  If  the  errors,  in  turn,  are  sums  of  a  con¬ 
stant  bias,  b,  plus  a  purely  random  error,  the  observations  may  be 
expressed  as  follows: 

R^(t)  =  Rj(t)  -l-b  +R  (t) 

where  R  (t)  is  the  random  error.  Subsection  1.3.2  shows  that,  with  a 
proper  choice  of  the  polynomial  degree,  the  curve  fitted  to  the  observa¬ 
tions  will  closely  approximate  the  sum  of  the  ideal  observations  and  the 
bias.  It  is  further  shown  in  subsection  1.3.  3  that  this  curve  will  not 
approximate  the  random  errors  very  well.  Thus,  the  technique  of 
polynomial  smoothing  separates  the  observations  as  follows:  the  curve 
fitted  to  the  observations  approximates  the  sum  of  the  ideal  observa¬ 
tion  and  the  bias;  the  residuals  or  differences  between  the  observations 
and  the  curve  approximate  the  random  errors. 

1.3.2  THEORETICAL  LIMITATION 

The  equations  for  the  polynomial  coefficients  are  linear,  thus  the 
estimated  coefficients  are  linear  functions  of  the  observations.  This 
implies  that  the  curve  of  best  fit  to  the  observations  is  given  by  the  sum 
of  the  curves  of  best  fit  to  the  ideal  observations,  bias,  and  random 
errors.  Therefore,  each  component  of  the  observations  may  be  analy¬ 
zed  separately. 

The  curve  of  best  fit  to  a  constant  bias  is  simply  that  bias  itself.  For 
the  ideal  observations, assume  that  a  polynomial  of  degree  "p  is  to  be 
fitted  to  range  observations  covering  a  time  span,  2T  .  Let  t^'  be 
the  midpoint  of  this  span.  The  ideal  range  in  a  truncated  Taylor 
Series  expansion  about  t^  is: 
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^  R,*'*  (t  )  (t-t  )"  R (f)  (t-t  )P*'' 

“i<''=Z  ^  ^  --■  Ip^y^--  . . 

r=0 

The  first  term  on  the  right-hand  side  can  be  fitted  exactly  by  a  poly¬ 
nomial  curve  of  degree  p.  The  second  term,  the  remainder,  involves 
the  (p+1) —  derivative  of  the  ideal  observation  evaluated  at  some  time, 

t',  between  t  and  t.  To  determine  how  closely  this  term  can  be  fitted 
o 

is  an  extremely  difficult  task,  however  its  magnitude  can  be  analyzed. 
The  magnitude  of  the  remainder  is  less  than; 


r(P+1) 

lpp+1 

I 

max 

(p+1): 


The  derivate  is  maximized  over  the  time  interval,  t-t^  <  T.  Thus, 
if  the  radius  of  convergence  of  the  Taylor  Series  is  greater  than  T, 
(Taylor  Series  converge  uniformly),  one  can  find  a  p  such  that  the 
remainder  terms  are  sufficiently  small.  In  practice,  p  can  be  chosen 
so  that  the  magnitude  of  the  remainder  term  is  a  small  fraction  (one- 
tenth  or  less)  of  the  expected  standard  deviation  of  the  observational 
errors.  Examples  of  the  magnitudes  of  these  standard  deviations  have 
been  computed  and  are  given  in  Appendixes  A  and  B. 

To  calculate  the  magnitudes  of  the  derivatives,  at  least  two  schemes 
are  possible.  They  could  be  computed  from  the  equations  of  ballistic 
motion.  In  the  LaGrangian  equations  of  motion,  the  second  derivatives 
of  range,  azimuth  and  elevation  are  expressed  in  terms  of  their  zeros 
and  first  derivatives.  An  estimate  of  the  zeros  and  first  derivatives 
can  be  obtained  from  the  record.  Differentiation  of  the  LaGrangian 
equations  the  proper  number  of  times  yields  the  desired  result.  This 
approach  rapidly  gets  out  of  hand,  however,  and  its  use  is  not  recom¬ 
mended. 
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The  second  scheme  is  more  simple  to  perform  but  is  subject  to 
numerical  computation  errors  from  an  assumed  orbit.  An  approximate 
orbit  from  the  record  is  determined.  From  this  orbit  the  radar  param¬ 
eters  at  equally  spaced  time  intervals  are  calculated.  The  derivatives 
of  these  radar  parameters  are  inferred.  Since  the  orbit  used  is  only 
approximate,  and  the  numerical  differentiation  introduces  errors,  the 
resulting  derivatives  should  be  used  with  a  safety  factor  on  the  maxi¬ 
mum  allowable  remainder. 

1. 3.  3  STATISTICAL  LIMITATIONS 

The  polynomial  coefficients  are  computed  on  the  basis  of  noisy  obser¬ 
vations  therefore  they  are  subjected  to  random  fluctuations.  For  the 
same  reason,  the  fitted  curve  is  also  subjected  to  fluctuations.  If  the 
fluctuations  in  the  curve  are  of  the  same  order  of  magnitude  as  the 
fluctuations  in  the  observations,  the  residuals  will  bear  little  resem¬ 
blance  to  the  observational  errors.  When  a  large  number  of  observa¬ 
tions  is  employed,  the  fluctuation  in  the  curve  becomes  quite  small. 

In  Appendixes  A  and  B  formulas  and  graphs  for  the  ratio  of  the  stand- 
dard  deviation  of  the  curve  to  the  deviation  of  the  observations  are 
presented  as  a  function  of  the  degree  of  the  polynomial,  number  of 
points,  and  time.  These  graphs  are  derived  on  the  basis  of  the  as¬ 
sumptions  that  the  observations  are  uncorrelated,  stationary,  and 
equally  spaced  in  time.  The  non- stationary  state  of  the  cartesian 
coordinates  does  not  make  them  strictly  applicable,  but  they  are 
nevertheless  approximately  correct. 

A  numerical  example  explains  the  reason  the  radar  coordinate  of  low¬ 
est  standard  deviation  gives  rise  to  distorted  residuals.  The  following 
table  shows  the  true  values  of  the  satellite  in  radar  and  cartesian  co¬ 
ordinates  and  the  radar  standard  deviations  at  a  particular  instant  of 
time: 


R 

E 

A  , 

■’R  "e  i  ‘'A 

X 

Y 

1 

800  1 

30° 

45°  ; 

2. 5  ,.0573°jo.l719°J 

490  j 

490  1 

Now  the  standard  deviations  in  X,  Y,  Z  may  be  found  as  follows: 

cos  E  cos  A)  ^  +  ( R  aj,  sin  E  cos  A)  ^  +  (  R  cos  E  sin  A^  ^ 
=  4.  744  =  (  2.  178)  ^ 

Similarly, 

=  {  2. 178)  ^  and  =(  1.429)^ 

Assume  next  that  the  standard  deviations  of  the  curves  fitted  to  the 
cartesian  coordinates  are  each  equal  to  one-twentieth  of  the  standard 
deviations  of  the  observations.  The  table  below  shows  these  deviations 
and  the  corresponding  deviations  in  the  curves  fitted  to  the  radar  coordi¬ 
nates: 


°Xc 

°Yc 

°Zc 

,  °Rc 

°Ec 

°Ac 

0.  109 

0.  109 

0.071 

i  0.101 

i 

0.00588° 

0.  00901° 

Then: 

‘^Rc/^^R  ”  ^Ec^^E  ■  0-103  t^Ac^^^A  ^  0.052 

Thus,  although  the  curves  were  fitted  equally  well  to  the  three  carte¬ 
sian  coordinates,  their  resulting  deviations  in  the  radar  parameter 
have  become  unequal  and  distorted. 
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1.3.4  DEGREE  OF  POLYNOMIAL 


A.  A  Priori 

A  priori  knowledge  of  the  degree  of  polynomial  required  to  fit  a  given 
set  of  observations  can  be  obtained  from  the  determination  of  polynomial 
fitting  to  the  true  orbital  path  for  the  given  length  of  observation  time. 
Although  the  true  trajectory  is  unknown,  a  family  of  trajectories  about 
the  ideal  value  can  be  investigated  to  determine  the  optimum  degree  of 
polynomial  required.  The  magnitude  of  the  errors  obtained  by  approxi¬ 
mating  the  ideal  trajectory  with  any  polynomial  is  indicated  in  Appendixes 
A  and  B. 

For  most  satellites  the  Millstone  Hill  radar  parameters  can  be  approxi¬ 
mated  by  a  fourth  order  polynomial  for  observation  times  up  to  5  4/2 
minutes.  Table  A-1  in  Appendix  A  indicates  the  time  spread  of  the  amount 
of  data  acceptable  for  various  deviations  in  the  fit  of  different  degree 
polynomials  to  the  ideal  path  of  satellite  motion. 

B.  A  Posteriori  (Orthogonal  Polynomials) 

Determining  the  polynomial  degree  may  also  be  done  a  posteriori  with 
the  use  of  orthogonal  polynomials.  The  polynomials  are  described  in 
Ref.  4o^d  only  certain  properties  will  be  mentioned  here. 

The  most  important  property  of  these  polynomials  is  that  they  may  be 
used  to  determine  the  polynomials  of  best  fit  of  several  different  degrees, 
with  a  minimum  of  duplication.  Closely  connected  with  this  property  is 
a  simple  statistical  test  (F-test)  which  may  be  made  on  the  difference 
between  the  variance  of  the  residuals  after  fitting  polynomials  of  succes¬ 
sive  degrees  to  the  data.  In  this  way,  an  a  posteriori  test  of  the  signifi¬ 
cance  of  the  p^^  order  polynomial  term  may  be  made. 


Use  of  orthogonal  polynomials  has  certain  computational  advantages 
over  the  usual  polynomials.  The  advantages  are  most  pronounced  when 
several  different  polynomial  degrees  are  to  be  tested.  These  advantages 
are  nullified,  if  an  a  priori  determination  of  the  degree  has  been  made 
or  if  the  observations  are  not  equally  spaced  in  time,  or  if  the  observa¬ 
tions  are  to  be  weighted  unequally. 

C.  Variate  Difference 

The  technique  of  applying  the  variate  difference  methods  to  the  observa¬ 
tions  is  considered  in  detail  in  Ref.  1.  Basically  it  determines  the  degree 
of  the  polynomial  fit  and  also  estimates  the  variance  of  the  random  com¬ 
ponent  by  considering  the  stability  of  the  forward  differences  of  the  param¬ 
eters.  It:assumes  that  the  radar  reports.are  composed  of  a  systematic 
and  a  random  component: , 

=^1  -t  "  n 

k=o 


where  p  is  the  order  of  the  time  polynomial  to  be  estimated.  Estimates 

of  the  variance  of  the  forward  differences  are  obtained.  This  variance 

til 

of  the  random  component  decreases  up  to  the  p  difference  and  then  the 
variance  should  stabilize  except  for  random  fluctuations.  These  random 
fluctuations  could  be  due  to  a  harmonic  periodicity  present  in  the  resi¬ 
duals. 

1,4  AUTOREGRESSIVE  SCHEME 

The  residuals  obtained  after  fitting  the  optimum  degree  of  polynomial 
still  exhibit  a  sinusoidal  tendency.  It  can  be  shown  that  the  fitting  of  the 
optimum  order  polynomial  plus  an  autoregressive  scheme  will  better 
describe  the  error  process.  In  addition,  this  method  gives  rise  to 
independence  between  the  observations  for  the  cases  considered  (Ref.  10), 
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1.4.1  DESCRIPTION  OF  AUTOREGRESSIVE  TECHNIQUE 

With  the  determination  of  the  order  p  of  the  polynomial  and  its  least 
square  coefficients,  the  residuals  can  be  fitted  by  an  autoregressive 
scheme.  The  autoregressive  scheme  considered  in  Ref.  10  is  of  the 
second  order: 


k=2 

k=l 

The  coefficients  b^  b^  2*  consistent  estimators  in  the  least 
square  sense.  In  addition  to  being  consistent  estimators,  these  least 
square  estimators  are  asymptotically  best  estimators.  The  solution 
to  the  autoregressive  difference  equations  has  a  harmonic  solution  which 
could  alter  the  behavior  of  residuals  obtained  from  polynomial  fitting 
only. 

The  solution  is  given  by: 

=  p*(Cj  cos  6t  +  C^  sin  dt) 

X,  =  Cj  B,*  t  s,‘ 

2 

depending  on  whether  the  roots  of  P  -b  P-b  ^  =  0  are  imaginary 

or  real. 

The  particular  solution  takes  the  form: 

j- 

Z  ^t-j+1 
j=0 

r 

i 
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so  that  initially  errors  are  independent.  Once  introduced  into  the  proc¬ 
ess,  they  forever  exert  their  influence.  For  this  reason,  it  was  ex¬ 
pected  that  residuals  would  be  uncorrelated,  the  dependence  tied  up  in 
the  initial  observations  and  remaining  uncorrelated  throughout  the  proc¬ 
ess,  and  therefore  tending  to  diagonalize  the  covariance  matrix.  An 
estimate  of  the  variance  of  the  autoregressive  scheme  is 


(Hb^)  (l-b^^  -  bj^) 


The  residuals  after  trend  and  autoregressive  fitting  were  considered 
elements  of  and  assumed  to  be  a  stationary  process  up  to  second 
order.  Whereas  only  the  second  degree  autoregressive  scheme  is 
indicated  in  Ref.  10,  the  degree  actually  required  can  be  estimated 
utilizing  similar  methods  to  those  used  to  determine  the  degree  of  the 
polynomial  via  a  variate  difference  method. 

1.5  EFFECT  or  REFERENCE  SYSTEM 

There  are  various  systems  of  coordinates  in  which  the  observation  can 
perhaps  be  smoothed  better  than  in  another  system.  The  systems  of 
radar  parameters  do  not  behave  like  low  order  polynomials  for  any 
length  of  time.  In  order  to  carry  out  polynomial  smoothing  on  these 
coordinates,  it  is  necessary  to  resort  to  polynomials  of  very  high  degree  - 
perhaps,  ten  or  twelve  or  more  for  a  typical  5-minute  pass.  It  is  not 
recommended  that  such  high  order  polynomials  be  used,  as  they  require 
the  solution  of  large  scale  systems  of  simultaneous  equations.  Such 
solutions  consume  much  time  and  storage  and  are  of  limited  accuracy. 
Another  objection  to  high  order  polynomials  is  the  loss  of  statistical 
stability  involved.  Furthermore  both  a  priori  and  a  posteriori  techni¬ 
ques  of  estimating  polynomial  degree  become  questionable. 
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Two  solutions  to  this  problem  are  known.  Either  a  sequence  of  functions 
may  be  found  (linear  and  nonlinear)  which  fit  the  observations  better  than 
polynomials  (autoregressive  scheme  perhaps)  or  the  observations  may  be 
transformed  to  a  form  in  which  they  are  easily  fitted  by  polynomials. 
Unless  differential  correction  is  carried  out  it  seems  quite  difficult  to 
measure  the  truncation  errors  associated  with  any  smoothing  other  than 
polynomial  fitting.  Furthermore,  most  smoothing  functions  will  give 
rise  to  nonlinear  systems  of  equations  which  will  involve  a  great  deal 
of  machine  computation. 

1.5.1  POLAR  VERSUS  TOPOCENTRIC  CARTESIAN  REFERENCE 
SYSTEM 

Although  the  radar  coordinates  range,  elevation  and  azimuth  do  not 
behave  like  low  order  polynomials,  the  cartesian  inertial  and  the  topo- 
centric  cartesian  coordinates  of  the  satellite  do.  Thus  the  transforma¬ 
tion  from  radar  to  topocentric  coordinates  produces  records  which  are 
very  amenable  to  polynomial  smoothing.  Appendix  A  shows  that,  for  a 
wide  variety  of  orbits,  the  truncation  associated  with  a  p  degree  poly¬ 
nomial  is  a  function  only  of  the  time  interval,  and  that  for  typical  radar 
data,  a  fourth  degree  polynomial  is  sufficient. 

The  weights  to  be  used  with  the  cartesian  coordinates  depend  upon  the 
standard  deviations  of  the  radar  coordinates.  Since  these  are  unknown, 
an  iterative  technique  is  employed.  This  technique  has  been  programmed 
for  the  IBM  7090  in  Fortran  language  and  performs  the  following  compu¬ 
tations. 

Initially  the  standard  deviation  in  range,  azimuth,  and  elevation,  and 
the  polynomial  degree  are  assumed.  The  observations  are  entered  and 
transformed  to  cartesian  coordinates.  Weights  for  the  cartesian  coordi¬ 
nates  are  found  on  the  basis  of  the  assumed  radar  deviations  and  the 
formulas  of  linear  error  propagation.  Then  normal  equations  for  the 
polynomial  coefficients  are  solved,  residuals  are  computed  and 


transformed  into  residuals  in  range,  azimuth,  and  elevation.  The  stan¬ 
dard  deviations  of  these  residuals  are  computed  and  compared  with  the 
assumed  values.  If  the  discrepancy  is  large,  the  procedure  is  repeated 
starting  with  the  computation  of  weights,  replacing  the  radar  deviations 
with  the  values  just  computed. 

This  program  has  been  applied  to  several  records.  In  all  cases,  con¬ 
vergence  has  been  achieved  after  two  or  three  iterations.  The  final 
residuals  have  not  deviated  significantly  from  the  Gaussian  (subsection 
1.6.  1).  Serial  correlations  have  been  small  for  those  coordinates  with 
relatively  large  standard  deviations,  but  have  been  significant  when  the 
deviations  are  small  -  usually  in  elevation.  This  need  not  contradict  the 
hypothesis  of  statistical  independence,  due  to  the  poor  statistical  stability 
in  these  residuals.  In  general,  it  is believedthat  the  residuals  in  those 
radar  coordinates  with  large  standard  deviations  are  close  approximations 
to  the  actual  random  errors,  while,  for  the  other  coordinates,  the  stan¬ 
dard  deviations  are  probably  reasonable,  but  the  residuals  themselves 
are  in  error.  An  example  of  the  statistical  limitations  and  distortion  of 
errors  is  indicated  in  subsection  1.  3.3. 

1.6  DATA  ANALYSIS 

1. 6.  1  RADAR  DATA  FROM  MILLSTONE  AND  MOORES  TOWN 

Samples  of  satellite  observations  from  two  different  sites.  Millstone 
and  the  FPS-49  radar  at  Moorestown,  were  analyzed  to  determine  the 
distribution  of  the  errors,  and  their  correlation.  The  data  supplied  by 
Moorestown  were  10-second  radar  parameters  for  various  satellites. 
These  satellites  included  1959  Epsilon  1  and  Iota  1,  I960  Zeta  1, 

Zeta  2,  and  Nu  2,  and  1961  Epsilon  1,  Lamba  1,  and  Omega  1.  The 
Millstone  data  were  6-second  semi-smoothed  radar  parameters  from 
I960  Epsilon  1  and  Iota  1.  Various  degrees  of  polynomials  were  fitted 
to  the  transformed  cartesian  coordinates.  The  statistics  of  the  resulting 
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residuals  of  the  radar  parameters  were  obtained.  From  the  satellites 
considered,  there  appeared  to  be  no  consistency  in  the  magnitude  of  the 
variances  of  the  individual  parameters.  A  typical  example  is  given  for 
I960  Epsilon  1  and  shows  the  effect  on  the  first  four  moments  of  fitting 
different  degree  polynomials  for  each  of  the  radar  parameters.  In 
general,  for  range,  azimuth  and  elevation  the  p  values  were  between  2 
and  3,  this  being  a  measure  of  the  normality  of  the  distribution.  The 
values  of  y  obtained  were  consistently  less  than  1,  i.  e. ,  the  residuals 
exhibited  a  tendency  toward  a  slight  skewness  to  the  right  for  the  bias. 

For  the  four  orders  of  polynomials  fitted  to  I960  Epsilon  1  the  covariance 
matrix  remained  relatively  unchanged.  It  can  also  be  seen  that  the 
covariance  was  not  a  diagonalized  matrix;  there  was  a  strong  correlation 
between  range  and  the  other  two  parameters.  A  usual  assumption  made 
about  radar  data  is  that  the  parameters  are  uncorrelated  to  facilitate 
easy  inversion  of  a  simple  diagonal  matrix.  Whereas  these  residuals 
with  a  polynomial  fitting  exhibited  a  marked  correlation,  when  the  radar 
parameters  are  fitted  with  a  polynomial  plus  an  autoregressive  scheme, 
as  indicated  in  Ref.  10,  the  covariance  matrix  obtained  can  be  assumed 
to  be  diagonalized.  Also,  the  standard  deviations  of  the  radar  param¬ 
eters  using  the  autoregressive  scheme  have  been  decreased.  The 
residuals,  after  fitting  the  autoregressive  scheme,  are  shown  in  Figures 
9  to  11  of  Ref.  10.  There  is  an  indication  that  a  periodic  trend  in  the 
residuals  is  still  present. 

The  serial  correlation  coefficients  of  the  radar  parameters  indicated 
that  the  residuals  obtained  were  independent  of  the  time  of  observation. 
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1.6.2  BAKER-NUNN  SATELLITE  TRACKING  CAMERA 


The  Baker-Nunn  tracking  camera  is  especially  designed  for  tracking 
artificial  earth  satellites  (Ref.  6).  The  camera  is  briefly  described 
as  a  three-axis  Super-Schmidt  f/ 1  camera,  with  a  focal  length  of  50  cm, 
and  a  field  of  view  5  by  30  degrees,  which  is  photographed  on  a  55-mm 
cinemascope  film  about  one  foot  long.  Faint  satellites  are  tracked  to 
obtain  measurable  images.  When  the  tracking  rate  of  the  camera  is  the 
same  as  the  rate  of  the  satellite,  the  satellite  appears  as  a  point  image 
on  the  film  against  a  star  trail  background.  The  position  of  the  satellite 
is  then  determined  by  measuring  the  centers  of  the  central  breaks  of 
reference  stars. 

The  number  of  observations  in  angular  position  that  the  Baker-Nunn 
camera  can  acquire  on  one  pass  of  a  satellite  depends  upon  the  satellite's 
altitude  and  brightness;  but  generally,  two  or  three  observations  are  taken. 
For  satellites  that  are  high  and  bright,  such  as  Iota  1,  it  is  possible  to 
obtain  ten  sample  observations  in  angular  position.  This  is  not  a  suffi¬ 
cient  number  of  observations  to  make  a  reliable  statistical  analysis. 

There  are  several  possible  schemes  for  obtaining  more  sample  observa¬ 
tions,  one  of  which  is  to  analyze  the  residues  of  the  observations  made 
on  many  passes  of  a  given  satellite  with  a  particular  Baker-Nunn  camera; 
but  each  pass  constitutes  a  different  experiment.  It  is  reasonable  that, 
with  each  experiment,  (in  effect,  measurements  under  different  geome¬ 
tries)  stationary  statistics  cannot  be  assumed  valid. 

A  second  scheme  that  might  be  suggested  is  an  ensemble  of  stations 
making  observations  on  a  given  satellite.  This  is  not  practically  feasi¬ 
ble.  In  addition,  in  analyzing  the  errors,  the  ensemble  averages  and 
the  time  averages  would  have  to  be  assumed  the  same.  In  either  scheme, 
the  process  is  neither  stationary  nor  ergodic. 
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Because  one  is  dealing  with  an  experiment  which  involves  so  few  sample 
points  and  is  devising  experiments  to  increase  the  sample  size, the  assump 
tion  of  stationarity  or  ergodicity  is  established.  However,  this  is  not 
valid.  The  approach  taken  by  the  Smithsonian  Institute  for  establishing 
accuracies  is  quite  therefore  reasonable  (Ref.  5). 

The  accuracies  established  for  the  Baker-Nunn  camera  were  determined 
as  follows:  Observations  from  stations  were  compared  with  respect  to 
a  very  accurately  known  orbit  and  the  residues  were  analyzed. 

Two  types  of  data  were  analyzed:  (1)  field  reduced  observations  and 
(2)  precisely  reduced  observations.  The  field  reduced  data  are  measure¬ 
ments  made  by  the  observers  at  the  stations  and  are  subjected  to  personal 
error.  The  field  reduced  positions  are  accurate  to  within  Z-l/Z  minutes 
of  arc  in  standard  deviation  and  one-tenth  second  of  time  in  standard 
deviation. 

Precisely  reduced  data  are  field  data  sent  via  airmail  to  the  Photo¬ 
reduction  Division  of  the  Smithsonian  Institute  in  Cambridge,  Massa¬ 
chusetts  and  analyzed  under  elaborate  laboratory  conditions.  Precisely 
reduced  data  have  an  accuracy  of  two  seconds  of  arc  and  0.  001  second 
of  time  in  standard  deviation. 

1.6.3  MULTI-SITE  INITIAL  WEIGHTING 

If  an  initial  knowledge  of  the  variance  of  the  data  can  be  estimated,  then 
it  can  be  utilized  for  quicker  convergence  in  the  Fortran  program  for 
the  polynomial  fitting  technique  as  described  in  subsections  1. 3.4  and 
1. 5. 1.  The  following  technique  has  the  advantage  of  obtaining  an  esti¬ 
mate  of  the  parameter  variance  of  one  site  if  the  parameter  statistics 
are  known  for  another  site.  The  assumption  made  is  that  the  variances 
are  unchanged  when  subjected  to  an  orthogonal  transformation.  To  carry 
out  this  weighting,  it  is  first  necessary  to  obtain  true  or  expected 
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parameter  values  for  azimuth,  elevation  and  range.  However,  at  any 
time,  for  any  satellite,  these  values  are  not  known.  In  Ref.  7,  equa¬ 
tions  are  given  for  the  determination  of  rectangular  coordinates  for  a 
satellite  in  a  geocentric  coordinate  system.  By  reworking  these  expres¬ 
sions,  the  following  procedure  is  obtained  for  prediction  of  the  satellite 
parameters  which  should  be  observed  at  one  radar  at  a  given  time,  t, 
when  the  parameters  are  known  for  some  other  radar  at  that  time,  t. 

The  known  inputs  are  the  positions  of  the  two  radar  sites  0^  6^, the 
geocentric  coordinates  of  each  radar  site,  and  the  radar  parameters 
from  both  sites  at  corresponding  times,  (S^^,  h^,  A^)  =  (slant  range, 
elevation,  azimuth).  (See  Figure  l-l.  ) 

The  subscript  merely  denotes  from  which  radar  site  these  parameters 
have  been  observed.  The  basic  equations  for  transferring  the  position 
coordinates  of  a  satellite,  as  observed  from  one  site  into  those  as 
observed  at  a  second  site,  are  given  as  follows; 


=  sin“^ 


(Rj  cos  hj) 


(Rj^  =  S(S+2Rj  sin 


"l  ■  2  ‘^^1 


The  geocentric  polar  coordinates  of  the  satellite  are  given  by 
r^  =  ^Rj^  +  S(S  +  2Rj  sin  hj)j 

J  ^  -  (cos  OTj  sin  +  sin  Oj  cos  cos  Aj)‘ 


and  also 


cos 


0=0,  -sin”^ 
s  1 


sin  OTj  sin 
cos  4 
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Figure  I'l.  Simplified  geometric  repreeentatioa  for  parameter 
prediction 


After  completing  the  transformation  the  geocentric  cartesian  coordi¬ 
nates  are  given  by 

> 

X  =  r  cos  6  cos  6 
s  s  s  s 

Y  =  r  cos  6  sin  6  , 

s  s  ^s  s  > 

Z  =  r  sin  ^ 
s  s  s 

* 

These  geocentric  cartesian  coordinates  can  now  be  converted,  as 
obtained  from  Site  1  data,  into  radar  parameters  as  would  be  observed 
at  Site  2. 


sin  -  cos  sin 

sin  a  cos 

s  2 

Thus,  by  considering  the  different  radar  sites  and  the  preceding  expres¬ 
sions,  values  from  one  radar  site  can  be  used  to  predict  corresponding 
values  for  another  site,  say  Millstone  Radar.  Also,  observed  values 
at  Millstone  will  be  used  to  predict  values  for  another  radar.  For  a 
given  time  interval  then,  four  sets  of  compatible  data  can  be  obtained. 
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For  a  given  time,  t,  an  observed  variance  S  ,  can  be  obtained  for  these 
data.  To  accomplish  this  the  following  information  is  considered. 

It  is  known  that  the  errors  between  the  observed  values  and  the  expected 
values  are  approximately  Gaussian  distributed,  N(0,  Op  ),  where  R  is 

Xx  “  p 

observed  value  and  p  is  the  expected  value. 

Thus , 


Thus, 

(R-p)  is  N(0,  a^) 

Z  Z 

For  this  reason,  with  S  known,  bounds  on  a  can  be  obtained  by  em¬ 
ploying  the  distribution  with  n-1  degrees  of  freedom  given  by 

y2  n 

X  =  - 2— 

a 

where  n  is  the  sample  size. 

2 

These  bounds  on  a  can  be  used  to  obtain  boundary  equations  (with  an 
attached  degree  of  confidence)  for  polynomial  smoothing  as  discussed  in 
subsection  1.3.  1 . 
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SECTION  2 


ORBITAL  COMPUTATIONS  AND  ERROR  BEHAVIOR 

Error  analysis  deals  with  the  type  of  errors  in  the  prediction  and  future 
ephemeris  implied  by  the  expected  errors  in  the  sensor  data.  This  con¬ 
cerns  the  mathematical  probability  that  the  vehicle  will  be  in  an  expected 
spacial  position  at  a  given  time. 

Error  behavior  can  be  studied  as  a  linear  or  a  nonlinear  problem.  The 
linear  analysis  is  based  on  assumptions  that  the  errors  in  the  prediction 
are  linearly  related  to  the  errors  in  the  input  data.  This  is  realistic  if 
the  standard  deviations  of  the  data  errors  are  extremely  small  and  pro¬ 
duce  prediction  functions  with  multivariate  Gaussian  distribution,  with 
their  covariance  matrices  related  by  simple  matrix  multiplication  in¬ 
volving  partial  derivatives  of  prediction  functions  with  respect  to  ob¬ 
servations.  The  way  the  error  analysis  is  performed  depends  upon  the 
distribution  of  statistics  assumed  and  the  dynamical  assumptions  applied 
to  the  arc  of  trajectory  during  the  period  of  observations.  The  first 
section  of  this  report  dealt  with  various  techniques  for  smoothing  these 
data.  This  section  discusses  the  error  behavior  of  the  osculating  ele¬ 
ments  and  finally  the  ephemeris. 

Statistical  assumptions  made  about  the  radar  parameters  are  that  they 
are  independent,  Gaussian,  and  can  be  expanded  in  a  Taylor  Series  and 
truncated  after  the  linear  terms  without  any  loss  in  accuracy.  Errors 
are  also  assumed  to  be  non-serially  correlated  and  time  invariant. 

While  several  of  these  assumptions  might  hold  for  some  radar  observa¬ 
tions  so  that  a  linear  error  analysis  could  be  used  for  the  osculating 


elements,  the  final  elements  would  certainly  not  be  independent  but 
would  exhibit  marked  correlation.  Some  results  obtained  in  the  pre¬ 
vious  section  tend  to  indicate  that  most,  if  not  all,  of  these  assumptions 
should  be  considered  open  to  suspicion.  The  only  assumption  that  seems 
to  be  allowed  is  that  the  time  of  observation  is  known  nearly  exact. 

2.  1  OSCULATING  ELEMENTS  DETERMINATION 

The  initial  point  of  the  study  of  error  behavior  is  the  conversion  of 
three  smoothed  vectors  into  one  smoothed  midpoint  velocity  vector  by 
a  method  attributed  to  Herrick-Gibbs.  This  technique  is  employed  by 
Space  Track  for  the  conversion  of  radar  data  (Ref.  7). 

An  alternative  method  is  to  obtain  the  velocity  vector  directly  from  the 
radar  data.  One  of  the  useful  by-products  of  the  topocentric  polynomial 
smoothing  is  a  measure  of  the  midpoint  three  components  of  velocity 
as  coefficients  of  the  first  order  terms  in  the  polynomial  and  the  accel¬ 
eration  as  the  second  order  terms.  The  statics  of  this  velocity  are  also 
easily  available.  Now  that  a  position  and  a  velocity  vector  are  both 
known  for  a  common  time^  the  osculating  elements  can  be  uniquely 
determined  (see  Appendix  D). 

The  smoothed  radar  data  exhibited  a  certain  amount  of  skewness  and 
non-normality  and  the  radar  parameters  were  definitely  not  independent. 
There  appears  to  be  no  serial  correlation  of  the  parameters  and  it  was 
not  ascertained  whether  the  covariances  were  time  invariant.  With  this 
knowledge  of  the  statistics  of  radar  smoothed  data  it  was  decided  to  in¬ 
vestigate  the  error  behavior  of  the  osculating  elements  (1)  on  the  basis 
of  linear  independent  errors  and  (2)  on  error  behavior  in  which  the  first 
four  moments  of  the  input  data  were  known  with  their  cocumulants. 

The  radar  parameters  were  ocpanded  in  a  Taylor  Series  about  their 
means  and  the  series  truncated  at  fourth  order  terms.  The  method  of 
propagation  of  errors  was  similar  to  that  of  Ref.  11.  It  quickly  became 
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obvious  that  this  was  an  unrewarding  tedious  task.  The  procedure  was 
reduced  to  only  the  first  two  moments  and  the  covariance.  The  neces¬ 
sary  partial  derivatives  and  formulas  for  the  error  prediction  of  the 
osculating  elements  are  indicated  in  Appendix  E. 

2.2  PREDICTION  DETERMINATION 

Error  behavior  of  the  final  prediction  was  analyzed  in  Appendix  E  by 
a  method  similar  to  that  employed  to  compute  the  osculating  element 
errors.  The  results  obtained  from  the  analysis  of  several  satellites 
indicated  that  the  osculating  elements  were  strongly  correlated  and 
were  not  Gaussianly  distributed  (see  subsection  2.  3).  Error  propaga¬ 
tion  for  the  first  four  moments  was  carried  out  to  evaluate  the  four 
moments  of  the  cartesian  coordinates  but  the  covariance  was  approxi¬ 
mated  to  second  moments  only.  Only  the  first  partial  differentiations 
were  carried  out  in  Appendix  E  due  to  their  lengthy  nature.  All  of  the 
partial  derivatives  have  a  periodic  tendency  which  is  non-decreasing 
with  time.  The  main  fault  of  lengthy  prediction  time  is  due  to  the  par¬ 
tial  derivatives  of  the  spatial  position  with  respect  to  the  semi-major 
axis  9r/8a. 

Although  this  derivative  was  also  periodic.- the  magnitude  of  this  period 
increased  with  time.  When  higher  order  derivatives  are  considered 
they  will  be  increasing  with  higher  orders  of  time.  Thus,  prediction 
error  will  increase  with  time.  It  might  be  predominant  compared  with 
the  errors  introduced  when  considering  a  simple  Keplerian  system, 
ignoring  the  effects  caused  by  earth's  oblateness,  drag  and  several 
other  second  order  effects  on  the  satellite. 

2.3  NUMERICAL  ERROR  BEHAVIOR  OF  OSCULATING  ELEMENTS 

Samples  of  Millstone  6-second  radar  data  for  three  different  satellites 
were  subjected  to  the  Herrick-Gibbs  technique  to  determine  a  series 
of  midpoint  velocity  vectors  throughout  the  observation  time  and  thence 
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the  resulting  osculating  elements.  Time  spacing  between  each  of  the 
three  position  vectors  selected  for  the  Herrick-Gibbs  technique  was 
varied  from  6  to  60  seconds.  This  increase  in  the  time  interval  con¬ 
siderably  reduced  the  variation  of  the  midpoint  velocity.  The  velocity 
vector  is, computed  from  three  weighted  position  vectors,  the  weights 
being 


The  time  spacing  is  such  that  G.  is  orders  of  magnitude  greater  than 
3  ^ 

H./y ^  which  is  the  component  that  considers  the  effect  of  the  Keplerian 
motion.  Since  both  G^  and  H.  are  functions  of  time,  this  magnitude  dif¬ 
ference  can  be  reduced  by  increasing  the  time  spacing.  The  longest 
time  spacing  indicated  from  the  Millstone  sample  data  is  approximately 

four  minutes  and  even  with  this  length  of  time,  the  G.  term  suppresses 

3  ^ 

the  Keplerian  motion  effect  of  hl^/y ^  . 

The  osculating  elements  were  then  fitted  with  unweighted  polynomials, 
the  order  of  which  was  predetermined  by  the  variate  difference  method 
and  the  statistics  of  the  residuals  of  these  osculating  elements  were 
derived.  The  complete  computation  procedure  from  radar  parameters 
-  osculating  elements  -  selection  of  degree  of  polynomial  and  statistics 
of  residuals  has  been  programmed  for  the  IBM  7090  in  Fortran  language. 

The  effect  of  increasing  the  time  spacing  of  the  data  had  the  effect  of 
decreasing  the  standard  deviation  of  the  osculating  elements,  shown  in 
Figure  Zrl  for  three  different  satellite  passes.  The  values  of  kurtosis 
varied  between  2  and  3  and  the  semi-major  axis  indicated  slight  skew¬ 
ness  to  the  right,  while  the  inclination  and  eccentricity  exhibited  the 
opposite  effect,  see  Figure  2-2,  2-3,  and  2-4.  Whereas,  the  cross 
correlation  for  the  radar  parameter  could  be  considered  negligible  for 
the  given  sample  size,  some  of  the  osculating  elements  have  a  strong 
correlation,  particularly  the  semi-major  axis,  eccentricity  and  inclina¬ 
tion  angle. 
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O  60  EPSILON  1  PASS  1701 
X  60  IOTA  1  PASS  167  SECTION  1 
^7  60  IOTA  1  PASS  167  SECTION  2 

Figure  2-2.  Kurtosis  and  skewness  behavior  of  eccentricity 
for  different  time  spacing. 


I 


1 


Page  33 
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X  60  IOTA  1  PASS  167  SECTION  1 
V  60  IOTA  1  PASS  167  SECTION  2 


6  12  18  24  30  36  42  48  54  60 


Figure  2-3,  Kurtosis  and  sKewness  behavior  ot  inclination 
lor  di  rent  time  spacing. 
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Figure  2-4.  Kurtosis  and  skewness  behavior  of  semi-major 
axis  for  different  time  spacing 


The  partial  correlation  coefficients  for  the  various  time  spacings 
appeared  to  be  relatively  stable  and  show  good  agreement  with  the 
simple  cross  correlation  for  both  the  radar  parameters  and  the  orbital 
elements.  The  double  sets  of  values  shown  for  1960  Epsilon  1  are  due 
to  the  correlation  between  two  parameters,  with  one  of  the  other  two 
parameters  fixed. 
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SECTION  3 


CONCLUSIONS 

This  study  has  shiown  that  the  differential  correction  smoothing  pro¬ 
cedure  is  a  solution  of  the  generalized  least  square  problem  (Ref.  2). 
The  residuals  obtained  by  fitting  the  observation  with  polynomial  plus 
an  autoregressive  scheme  (Ref.  10)  are  uncorrelated  and  the  standard 
deviations  of  the  radar  parameters  with  the  autoregressive  scheme  are 
less  than  those  obtained  with  polynomial  fitting  alone. 

In  the  process  of  analyzing  the  behavior  of  the  orbital  elements  it  was 
noted  that  the  coefficient  of  skewness  exhibits  a  degree  of  consistency 
for  different  satellites  and  sensors.  Furthermore  it  has  been  possible 
to  derive  the  variance  of  any  orbital  elements,  as  functions  of  the  vari¬ 
ance  of  the  observations,  the  time  of  observation  and  description  of  the 
motion  of  the  satellite.  It  has  been  shown  that  the  variance  of  any 
orbital  element  is  a  minimum  if  the  sum  of  the  variances  of  the  obser¬ 
vations  at  a  given  time  is  a  minimum  for  a  specific  computation  pro¬ 
cedure  of  the  orbital  element  from  observations.  This  indicates  that 
minimum  variance  orbital  estimations  are  only  relative  minimums. 

The  behavior  of  the  variance  and  skewness  suggests  that  the  computa¬ 
tional  procedure  introduced  a  bias.  The  latter,  in  particular,  cannot 
be  traced  to  a  round-off  error. 

The  problem  of  getting  accurate  position  prediction  of  artificial  satel¬ 
lites  is  one  of  interplay  between  the  quality  of  observations  and  the 
analytical  description  of  the  motion  of  the  satellite.  How  to  trade  these 
factors  depends  upon  the  time  allowed  for  deriving  the  orbit  from  obser¬ 
vations. 
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The  analytical  expression  showing  this  interplay  has  been  derived  and 
represents  the  natural  starting  point  for  (1)  an  accurate  prediction 
scheme  and  (2)  updating  or  intermixing  criteria. 
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SYMBOLS  USED 


p(t^) 

A(t.) 

E(t.) 

<t> 

X 

t. 

(  X,  y,  z) 


r(  t.) 


R 


r(  t.) 


d.,  G.,  H. 

Ill 

2 

s 


P 


U,  V,  w 
a 


e 

V 


a 

CJ 


N,  M 
G 


Range  at  time  L 

Azimuth  at  time  t. 

1 

Elevation  at  time  t. 

1 

Site  latitude 
Site  longitude 
Time 

A  geocentric  coordinate  represented  vector 
(x,  y,  z)  in  topocentric  coordinates 

A  vector  from  earth's  center  to  a  point  on  orbit  (at  L) 
Earth's  radius 
Time  increments 
Rate  of  change  of  r(t.) 

Defined  in  text,  functions  of 
Defined  in  text 
Defined  in  text 
Defined  in  text  (vectors) 

Semi-latus  rectum  of  orbit 
Eccentricity  of  orbit 
True  anomaly  of  orbit 
A  vector,  defined  in  text 
The  angle  of  perigee 
Vectors  defined  in  text 
Gravitational  constant 
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E 

T 

Q 

I 


I,J,K 

E(x) 

P 

P 

P 

a 

xy 


Px 


xyz 


xyzo) 

A,  B 


C 

xy 

C 


xyz 

C 


xyzw 


Eccentric  anomaly  of  orbit 
A  constant  of  the  orbit 
The  angle  of  the  ascending  node 
The  inclination  of  the  orbit 
Unit  vectors,  defined  in  text 
The  expectation  of  X 

The  "ith"  moment  of  P  about  its  mean,  i=2,  3,  4 
The  mean  of  P 
E(  x-x)  (y-y) 

E(  x-x)(y-y){z-z)/a^a^a ^ 

E(x-x)(y-y)(z-z)(w-w)  / 

Matrices,  3x3 
An  element  of  A  B 
E(x-x)(y-y) 

E(x-x){y-y)(z-z) 

E(x-x)(y-y)(z-z)(w-w) 
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APPENDIX  A 


POLAR  VERSUS  CARTESIAN  SMOOTHING 

Figure  A-1  is  a  rough  graph  of  the  fourth  degree  curve  of  best  fit  to  ten 
minutes  of  Millstone  azimuth  data.  It  is  important  to  note  that  although 
the  observations  cluster  around  a  smooth  curve,  this  curve  differs  con¬ 
siderably  from  the  fitting  curve.  This  is  because  radar  azimuth  and  ele¬ 
vation,  and  to  a  lesser  extent,  range,  do  not  behave  like  low  order  poly¬ 
nomials  for  any  length  of  time.  To  adequately  represent  the  observations 
by  a  fourth  degree  polynomial,  it  is  necessary  to  restrict  the  data  to  a 
much  smaller  interval,  perhaps,  one  to  two  minutes.  At  the  Millstone 
repetition  rate,  this  allows  only  ten  to  twenty  pulses  to  be  smoothed. 

As  it  will  be  seen,  a.  larger  number  of  pulses  is  needed  if  the  resulting 
observational  errors  are  to  have  any  statistical  meaning. 

Although  the  radar  coordinates  do  not  behave  like  low  order  polynomials, 
the  cartesian  inertial  coordinates  of  the  satellite  do.  The  amount  by 
which  these  coordinates  deviate  from  a  polynomial  of  degree  p,  during 
a  time  interval  of  length  2T  is  certainly  less  than  the  following  expres¬ 
sion; 


(Dev)  < 
P 


xP+1 


r(P+l) 


max 


(A-1) 


where  is  the  p+1*^  derivative  of  the  geocentric  position  vector, 

the  double  bars  indicate  its  magnitude,  and  "max"  refers  to  the  maxi¬ 
mum  value  of  this  magnitude  over  the  time  interval.  Equation  (A-1) 
is  an  application  of  Taylor's  theorem. 
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In  Table  A-1,  numerical  bounds  are  found  on  the  quantities  max 

by  noting  that  (1)  these  quantities  may  be  calculated  without  signi¬ 
ficant  error  assuming  a  Keplerian  ellipse,  and  (2)  on  this  ellipse  the 
derivatives  are  maximized  at  perigee.  The  region  of  a-e  space  indicated 
in  the  third  column  includes  most  reasonable  and  stable  satellite  orbits. 
All  formulas  and  numbers  are  in  canonical  units. 

Table  A-1 


Magnitude  of  first  few  position  derivatives  is  evaluated  at  perigee  and 
numerical  bounds  on  these  magnitudes. 


TT 


4  (143e)r/5’  2.0 

5  VlTc  (l+9e)  r  ^  4.  8 

TT 

6  (  l+24e-f  27e^)  r  °  13.  0 

TT 
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The  second  or  third  column  of  Table  A-1  may  be  used  with  Equation  (A- 1) 

depending  upon  whether  the  orbital  elements  are  approximately 

known.  Combining  the  third  column  with  Equation  (A-1),  Table  A- 2 

may  be  obtained.  Table  A- 2  expresses  the  maximum  values  of  the  in- 

til 

tervals  over  which  the  inertial  coordinates  deviate  from  a  p  degree 
polynomial  by  the  given  amounts. 


Table  A-  2 


2T  versus  p  and  (Dev) 


P 


p 

(Dev)  <  1 . 
P 

1  km 

(Dev)  <  0, 
P 

1  km 

(Dev)  <  0. 
P 

'  1 
1 

01  km 

0 

0.205 

sec 

0.0205 

sec 

0. 00205 

sec 

1 

29.9 

sec 

9.  56 

sec 

2.99 

sec 

2 

157. 00 

sec 

72.  90 

sec 

33.8 

sec 

3 

334.00 

sec 

188.00 

sec 

106. 00 

sec 

4 

529.00 

sec 

334.00 

sec 

210. 00 

sec 

5 

726.00 

sec 

499.00 

sec 

337. 00 

sec 

Table  A- 2  shows  that,  if  a  0.  1-km  deviation  in  the  fit  can  be  accepted,  a 
fourth  degree  polynomial  may  be  applied  to  5-1/2  minutes  of  data. 

Thus,  the  following  procedure  has  been  developed  for  the  orbit  removal 
from  Millstone  radar. 

(1)  From  each  pass,  select  the  middle  seven  minutes  of  data 
(Dev^  (7  min)  <0.35  km). 

(2)  Transform  the  ranges,  azimuths,  and  elevations  selected  in 
Step  (1)  to  Cartesian  inertial  coordinates  by  the  well-known 
transformatio 
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(3)  Solve  the  three  sets  of  likelihood  equations  for  the  selected 
records  of  x,  y,  z. 

(4)  Use  the  solutions  to  the  likelihood  equations  to  evaluate  the 
fourth  degree  polynomials  for  x,  y,  z  at  all  the  tj^  points 
selected  in  Step  (1). 

(  5)  Transform  these  polynomial  points  back  to  range,  azimuth 
and  elevation  and  subtract  them  from  the  corresponding 
observations. 

One  simplification  of  this  scheme  is  possible.  Step  2  may  be  replaced 
by  a  transformation  to  the  cartesian  non-inertial  coordinate  system; 

X  =  R  sin  E 

y  =  R  cos  E  sin  A  (A-2) 

z  =  R  cos  E  cos  A 

This  procedure  is  simpler  than  the  inertial  transformation.  The  justi¬ 
fication  of  this  simplification  is  given  by  noting  the  relationship  between 
the  usual  inertial  system  and  the  system  defined  by  equations  (A-2)  as 
shown  in  Figure  A-2.  The  inertial  system  has  its  origin  at  the  earth's 
center  with  a  z-axis  along  the  earth's  north  pole,  an  x-axis  pointed  to 
the  first  point  of  Aries  (y),  and  a  y-axis  completing  the  right-hand  sys¬ 
tem.  The  topocentric  system  has  its  origin  at  the  position  of  the  radar; 
its  x-axis  points  directly  away  from  the  earth's  center;  its  z-axis  points 
in  the  direction  of  true  north.  Let  the  vector  position  of  the  satellite 
in  the  two  systems  be  r,  and  p,  respectively,  then  the  transformation 
is  obtained: 

p  =  -Cj  +  ABr  (  A-3) 

The  matrices  A  and  B  and  the  vector  are  given  by: 
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Figure  A-2.  Topocentric  and  inertial  coordinate  aysteme 
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/  cosa 
=  (  -sina 
\  0 


0 

0 

1 


sina  \ 
cosa  I 
0  / 


B  = 


sinP 

cos^ 

0 


(A-4) 

In  Equation  (A-4),  a  is  the  declination  of  the  observer,  the  P  is  his 
right  ascension  at  any  time.  By  differentiating  Equation  (A- 2): 


=  ABr^^^  +  5ABr^^^  +  lOABr^^^  +  lOAB^^^r 


+  5AB^‘^^  r  +  AB^^^  r 


(A- 5) 


Noting  that  B^^^  is  equal  to  times  a  suborthogonal  matrix  (a  rotation 
plus  a  projection)  then: 


if  5P 


<  4.8  +  5(0.0587)  (2.0)  +  10(0.0587)^  (1.0)  +  ... 


(A-6) 


<5.4 


This  is  a  somewhat  conservative  bound,  probably  4.8  is  more  realistic. 

In  any  case,  a  fourth  degree  polynomial  fits  the  data  equally  well,  regard¬ 
less  of  whether  the  data  are  expressed  in  inertial  or  topocentric  coordi¬ 
nates. 
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APPENDIX  B 


DISTORTION  OF  THE  ERRORS 


Since  the  polynomial  coefficients  estimated  from  a  record  are  random 
variables,  the  calculated  values  of  the  polynomials  at  each  t^  are  also 
random  variables.  If  these  values  have  standard  deviations  as  large 
as  the  deviations  of  the  observations,  the  actual  observational  errors 
will  be  masked.  The  deviations  of  the  calculated  values  are  calculated 
here  under  two  assumptions: 


(1)  or  X/4.j\  =0'^  611,'  That  is,  the  data  have  no  serial 


2 

x(  tn)  *(tj)  ■  “  X  °  jk' 
correlation,  and  the  same  coordinate  observed  at  different 
times  has  the  same  variance. 


(2)  tj^^j  =  tj^  +  T.  That  is,  the  observations  are  equally  spaced 
in  time. 


Assumption  (1)  is  not  valid,  for  even  if  it  were  true  of  range,  azimuth, 
and  elevation,  it  would  not  be  true  that  the  cartesian  coordinates  would 
have  constant  variances.  Assumption  (2)  is  approximately  satisfied 
for  any  radar  -  in  particular  for  Millstone,  T  =  6  seconds. 

Given  assumption  (1),  Equation  (B-1)  for  the  covariance  matrix  of  the 
coefficients  is  valid  (with  replaced  by  <7^,  or  a^).  Given  assump* 
tion  (2),  the  sums  in  the  matrix  may  be  calculated  in  closed  form,  and 
the  matrix  analytically  inverted.  The  results  of  this  inversion  are 
shown  in  Table  B>1.  The  left-hand  column  is  the  degree  of  the  polyno¬ 
mial.  The  matrix  on  the  right  hand  has  o  .  ,  A.  ,  in  its  i^  row  and 

rt,  •'n-1 

j  column.  The  quantity  T,  as  defined  previously,  is  one-half  the  time 
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interval,  and  the  polynomial  coefficients  are  evaluated  at  the  time  corres¬ 
ponding  to  the  midpoint  of  the  interval. 


4 


2 

A. 


^0^1 


AA*^  !  (Ta  A 

i  -^0-^1  -^1 


A  A 


\ 


A  A 


•  •  •  ^ 


^2^4 


y 


=  <Tt 


^4) 


=4 


J 


These  covariance  coefficients  are  obtained  from  the  maximum  likelihood 
equation  of  observation  that  are  random  variables. 

Table  B-1 

Approximate  covariance  matrix  of  midpoint  polynomial  coefficients  ver¬ 
sus  degree  of  polynomial,  as  a  function  of  number  of  points,  N;  half¬ 
time  interval,  T;  and  variance  of  observations,  a^. 


pa™ 

,  AA^ 

0 

1 

N 

1 

/  N 

12^0 
jy  a  Q  3/.J.2 

V  > 

1 

r- 

3  0  -5/T^ 

0  4/T^  0 

-5/T^  0  15/T^ 

k  > 
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Table  B-1  (Continued) 
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0  -210/T^ 

0 

240/T^  0 

-210/T^ 

0  1764/T^ 

0 

-336/T^  0 

189/T^ 

0  -1890/T 

0  189/T^ 

-336/T^  0 
0  -1890/T^ 

560/ 0 
0  2205/T® 


The  data  in  Table  B-1  are  inaccurate  for  N  <  '^  20.  For  N  >  20,  they 
are  in  excellent  agreement  with  the  true  matrices. 

From  these  data,  the  variances  of  the  polynomial  points  may  be  deter¬ 
mined.  For  example,  for  the  point: 


(t)  •  Ap  +  A,  t  +  t  Aj  t  A^  t* 


(B-2) 


The  variance  is  obtained  by  the  formula: 


(B-3) 
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as: 


%(t)  = 


3.  52  -  14.06  +  114.84 

f  171.95  (ijt)® 


-  251.  56  (y)^ 


{B-4) 


For  polynomials  of  lower  degree,  expressions  similar  to  Equation  (B-4) 
are  found.  These  expressions  are  listed  in  Table  B-2.  In  Figure  B-1 
the  square  roots  of  these  expressions  and  the  ratio  of  the  corresponding 
standard  deviations  are  plotted  for  several  degrees  and  several  values 
of  N. 


Table  B-2 

Ratio  of  variance  of  polynomial  point  to  variance  of  observation, 

V  various  degrees  of  polynomial,  p,  as  a  function  of  N 

V ' 

and  of  the  ratio  of  time  from  midpoint  to  half  interval,  K  -  t/T. 


p 

/„2 
®Xp(t)' ^x 

0 

a^/N 

1 

oVn  ( 1+3K^) 

2 

oVn  (  2.  25-4K^+U.  25K^) 

3 

oVn  (  2.  25+1 1 . 25K^-41.  25K^+43.  75K^) 

4 

oVn  (3.516-14.064K^+114. 844K^-251.  562K^+71. 953K®) 

Figure  B-1  indicates  that  with  Ns 70  and  fourth  degree  smoothing,  up  to 
80  percent  of  the  xiomputed  points  have  a  standard  deviation  less  than 
one-quarter  of  the  observation  deviation;  furthermore,  at  the  extreme 


I 


54 
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Figure  B-1.  Ratio  of  standard  deviation  of  estimated  polynomial 
1  point  to  standard  deviation  of  observations,  versus 

!  ratio  of  time  from  midpoint  to  half  interval,  for 

second,  third,  and  fourth  degree  smoothing. 
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edges  of  the  observation  interval,  this  ratio  is  about  three-fifths.  How¬ 
ever,  since  assumption  (1)  of  the  previous  section  -  that  the  observa¬ 
tion  deviations  in  each  coordinate  are  constant  -  is  not  true,  the  ratio 
may  be  somewhat  higher  than  this  value.  For  this  reason,  it  is  con¬ 
sidered  advisable  to  chop  off  the  errors  at  about  80  percent  of  the  way 
to  the  edges. 

A  smoothing  procedure  has  thus  been  established  to  be  applied  to  Mill¬ 
stone  data  covering  seven  or  more  minutes: 

(1)  Select  the  middle  seven  minutes  of  data. 

(2)  Transform  range,  azimuth  and  elevation  to  cartesian  coordi¬ 
nates  by  Equations  (A- 2). 

(3)  Find  the  coefficients  of  the  fourth  degree  polynomials  of  closest 
fit  to  the  three  records  obtained  in  Step  (2)  by  solving  the  likeli¬ 
hood  equations  assuming  observations  are  random  variables. 

(4)  Evaluate  the  three  polynomials  found  in  Step  (3)  at  each  of  the 
71  time  points  selected  in  Step  (1). 

(5)  Transform  the  polynomial  points  found  in  Step  (4)  to  ranges, 
azimuths  and  elevations  by  inverting  Equations  (A- 2)  and  sub¬ 
tract  these  values  from  the  observed  values  of  range,  azimuth, 
and  elevation. 

(6)  Discard  the  first  and  last  seven  sets  of  differences  leaving  57 
"range  errors,"  57  "azimuth  errors,"  and  57  "elevation 
errors. " 

This  procedure  has  been  devised  for  Millstone  data,  which  are  availa¬ 
ble  at  6-second  intervals,  and  which  have  standard  deviations  of  about 
one  to  two  kilometers.  For  different  radars,  a  different  procedure  is 
necessary,  but  the  formulas,  tables,  and  graphs  presented  here  should 
be  su^icient  to  derive  such  a  procedure. 
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APPENDIX  C 


CONFIDENCE  BOUNDS  FOR  SAMPLE  CORRELATION 
COEFFICIENTS  AS  A  TEST  FOR  RANDOMNESS 


In  order  to  determine  a  test  of  randomness  for  a  set  of  samples,  an 
initial  hypothesis  p  =  0  is  taken.  On  the  basis  of  the  serial  correla¬ 
tions  of  the  sample,  a  procedure  for  testing  randomness  by  employing 
the  distribution  function  of  the  correlation  coefficients  is  established. 

A  90-percent  confidence  curve  on  the  correlation  coefficient  is  shown  in 
Figure  C-1 . 

Given  a  sample  of  size  m,  it  is  desired  to  test  the  sample  for  random¬ 
ness.  One  approach  to  this  problem  rests  upon  an  examination  of  the 
serial  correlations.  The  serial  correlation  coefficient  is  given  by 


COV  (x.,  X.  ) 

— J.— 


var  X., 
J 


var  X 


•j+k 


'V72 


k  =  0,  1,  2  . . . 


(C-1) 


If  the  elements  x^,  are  separated  into  two  groups,  then  the  ele¬ 

ments  X.  and  X.,,  can  be  considered  as  elements  from  separate  sam- 
pies,  both  samples  being  of  the  same  size,  n. 


It  is  now  assumed  that  the  population  has  a  bivariate  normal  distribu¬ 
tion.  That  is,  the  samples  x.  and  x.  ,  can  be  considered  as  coming 

3  jtk 

from  a  bivariate  normal  population.  The  distribution  function  of  r  as 
given  in  Ref.  8  is: 
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vermas  n 


Figure  C-1.  90  percent  confidence  curve  for  serial  correlation 

coefficients  for  various  sample  -  sizes. 


(C-2) 


q(r) 


1- 

n-l/2 

(n-)/2 

]-r 

'n-zr 

L  2  J 

X  ^  [(2pr)^/^'J  ‘p  [(n-l+k)/2|,  -1  <  r  <  1,  p  =  Er 

k=0 


for  large  n,  Er  =  p,  a  = 


-  U-pV 


n 


Tables  of  the  distribution  q(r)  can  be  found  in  Ref.  9. 


(C.3) 


A  test  for  randomness  under  the  initial  hypothesis  is  now  considered. 
Tor  p  =  0,  q(r)  reduces  to 


q(r) 


n-4 

2 


(C-4) 


Solving  this  expression  for  r, 


r 


VT  qCr)!-*! 

r  i=r) 


n-2 

2 


2 

n-4 


(C-5) 


For  a  given  q(  r),  n,  bounds  can  be  established  on  r,  <  r  <  r^,  in 
order  to  test  p  =  0  with  any  level  of  confidence.  The  bounds  may 
be  rewritten  as  /r/  <  r^  since  r  can  be  seen  to  be  symmetric  about 
zero. 
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To  use  the  tables  in  Ref.  9i  values  of  p  and  n  are  needed.  For  a  given 

probability,  a  value  of  r  can  be  determined.  One  such  set  of  values  is 

included  in  this  report.  These  values  are  for  a  Two-Tailed  Test  of  H  : 

o 

p  =  0  with  a  probability  of  0. 90.  The  curve  which  follows  the  table  is 
merely  a  graphical  representation  of  the  results  with  only  the  upper 
curve  recorded  since  symmetry  exists  about  zero. 

For  each  serial  correlation  coefficient,  the  number  of  pairs  of  observa¬ 
tions  is  plotted  against  r  on  the  graph.  If  |r|  falls  under  the  curve,  it 
may  be  saidthat  the  sample  is  random  with  a  probability  of  0.  90. 


60 


I 


CRITICAL  VALUES  OF  fi  FOR  90  PERCENT 
CONFIDENCE  BOUNDS 


D,  P  =  0.05 


IT,  P  =  0.95 


-0.975267 

-0.900 

-0.805676 

-0.730231 

-0.670540 

-0.622746 

-0. 583502 

-0. 544231 

-0.522944 

-0.497620 

-0.477867 

-0.458413 

-0.441914 

-0.427845 

-0.413950 

-0.400026 

-0.390138 

-0.380474 

-0.370914 

-0.361411 

-0.351838 

-0. 344748 

-0.338367 

-0.238043 

-0.170060 

-0.117231 

-0.083144 


0.965267 

0.900 

0.805676 

0. 730231 

0.670540 

0.622746 

0.583502 

0. 544231 

0.522944 

0.497620 

0.477867 

0.458413 

0.441914 

0.427845 

0.413950 

0.400026 

0.390138 

0.380474 

0.370914 

0.361411 

0.351838 

0.344748 

0.338367 

0.238043 

0.170060 

0.117231 

0.083144 


APPENDIX  D 


DERIVATION  OF  THE  ORBITAL  ELEMENTS 

This  appendix  summarizes  the  equations  needed  to  obtain  the  orbital 
elements  from  observation  data.  The  Herrick-Gibbs  method  is  used  in 
determining  the  position  vector  and  velocity  vector  of  an  observed  satel¬ 
lite. 

D.  1  DETERMINATION  OF  INERTIAL  COORDINATES 

In  general,  the  observations  of  a  satellite  are  made  with  reference  to  a 
given  site  of  known  latitude  and  longitude.  To  derive  the  orbit  of  this 
satellite  the  observed  positions  of  the  satellites  are  referred  to  as  a 
geocentric  inertial  reference. 

Let  the  earth's  center  be  the  origin  of  the  geocentric  inertial  reference. 
The  axis  of  the  earth  represents  one  of  the  axis,  Z,  of  the  reference. 

One  of  the  axfes,  X,  is  in  the  equatorial  plane  and  directed  toward  the 
first  point  in  Aries  or  vernal  equinox;  the  other,  Y,  is  also  in  the  equa¬ 
torial  plane  and  orthogonal  to  X.  Let  X',  Y',  Z'  denote  a  topocentric 
reference  system  whose  origin  is  the  site.  The  X'  £ucis  coincides  with 
the  vertical  through  the  site,  Z*  is  the  horizontal  plane  about  the  site  and 
directed  toward  the  north,  and  Y'  forms  an  orthogonal  system. 

Let  0  be  the  angular  rate  of  the  earth's  rotation  in  radians  per  second 
and  let  the  earth  be  spherical  of  radius  R.  If  X',  Y',  Z'  is  a  point  in 
the  topocentric  reference  and  X,  Y,  Z  is  a  point  in  the  inertial  reference 
then: 
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”1 


7  V  • 

where  ^  is  the  site  latitude,  X  is  the  site  longitude  and  X  =  X  +  ^  t. 

If  the  observations  of  a  satellite  about  the  site  S  are  range,  p,  elevation, 
E,  and  azimuth,  A, 


P  X’  (t}\  ^  sin  E 

P  y'  (t)  }-(  p  cos  E  sin  A 
z’  (t)y/  \p  cos  E  cos  Ay 


then  the  position  of  the  satellite  in  the  inertial  coordinate  is  represented 
by 


(cos  ^ 

sin  ^ 


0  -sin  A  iiitj)  sin  E(tj) -t-R 
1  0  I  p(tj)  cos  E(tj)  sin  A(tj) 

0  cos  \p(tj)  cos  E(tj)  cos  A(tj)/ 


where  t^  denotes  the  time  of  the  observation 

D.2  DETERMINATION  OF  POSITION  AND  VELOCITY  VECTOR 

From  observations  through  vector  positions,  r  (tj),  r  (t2),  £  (tj)  are 
determined  using  the  transformations  in  subsection  D.  1.  These  posi¬ 
tions  are  determined  at  approximately  equal  intervals  of  time.  With 
these  vectors  r  (t^)  and£(t2)  are  derived  using  the  Herrick-Gibbs 
method. 
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r  ( t2)  =  I  (  ^2) 


r  (t2)  = 

-dj  £  ( tj)  +  d2  r  (t2)  +  <13  £  (^3) 

where 

d.  =  Gi  f  H,  [r  (tj)]  -5 

and 

Hj  =  T2/I2.  H2  =  (T2  -  Tj)/12. 

II 

X 

[VjT2(Tj  tTj)]  O^MT^-rf) 

'^2' ^1*^2 

=  Tf[^TiT2(TjtT2)]  ■* 

where 

T,  andTj^tj-tj 

D.  3  DETERMl  NATION  OF  ORBITAL  ELEMENTS 

With  the  values  of  r  (t2)  and  r  (t2)  the  orbital  elements  are  determined 
as  follows: 

2 

r  *  £  *  £ 

•  2 

s  =  £  •  £ 
r  =  r •  r/r 
4  £=r]r|-^ 

r  r  -  r  r 


V  = 


W  *  U  X  V 


1  -  2 

a  ■  "FT 


2  1 
e  =  1 


2, -2  *2. 

r  (s  -  r  ) 


r  1-2  ,  2  ^ 

e  8in  v  =  7T  Vr  (s  -  r) 


e  cos  V  =  - i - jm -  -1 


V  =  tan 


■  ^  fg  8in  V  1 
[e  cos  vj 


U  e  cos  V  -  V  e  sin  v 

X  X 

a  =  U  e  cos  v  -  V  e  sin  v 

-  y  y 

U  e  cos  V  -  V  e  sin  v 
z  z 


a  =(-aW+aW)//l-W^ 
xn  X  y  y  X  '  V  z 


=  (-a  WW  -  aWW)/  /  1  -  +  a  /  1  -W^ 

'xxz  yyz'v  z  zV  z 


«  =  tan”  — ^  ,  u  =  «  f  V 

xn 
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U  cos  u  -  V  sin  u 

X  X 

N  =  U  cos  u  -  V  sin  u 

-  y  y 

U  cos  u  -  V  sin  u 
z  z 


-  U  sin  u  +  V  cos  u 

X  X 

M  =  U  sin  u  +  V  cos  u 

y  X 

U  sin  u  +  V  cos  u 
z  z 
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( 


I 

] 

1 

1 

1 

] 


Thus  the  six  orbital  elements  are  defined, 
[a.  e,  1.  <ji,  fi.  t] 


i 

I 

] 

) 

I 
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(1,  0.  0) 


=  (0.  0.  1) 


namely 


APPENDIX  E 


THE  MOMENTS  OF  PREDICTED  POSITION  ERROR 
FROM  MOMENTS  OF  ORBITAL  ELEMENT  ERRORS 

This  appendix  presents  the  nnethod  for  obtaining  the  statistics  of  the 
error  in  the  predicted  position  of  a  satellite  as  a  function  of  the  statis¬ 
tics  of  the  error  in  orbital  elements. 

E.  1  CALCULATION  OF  INERTIAL  GEOCENTRIC  COORDINATES 

FRTSKTQgBITAE  PARAKIETER - 

E.  1 . 1  STATEMENT  OF  THE  PROBLEM 

The  determination  of  ephemeris  of  particular  satellites  from  radar 
data  which  are  subject  to  various  errors  will  give  predicted  parameters 
subject  to  errors.  The  inputs  are  the  orbital  elements,  which  have  as¬ 
sumed  moments  and  comoments.  The  geocentric  and  topocentric  iner¬ 
tial  coordinates  are  then  derived  with  their  respective  moment  distribu¬ 
tions.  The  initial  orbital  information  was  considered  to  be  non-indepen¬ 
dent  data,  and  the  analysis  utilized  the  first  four  moments  to  derive 
formulas  for  the  corresponding  moments  of  the  geocentric  inertial  coor¬ 
dinates.  These  formulas  indicate  the  procedure  for  the  computation 
of  such  terms,  but  due  to  their  length  they  were  approximated  to  include 

7 

only  powers  of  order  <  a  as  they  appear  in  Equations  E.  2. 12,  E.  2. 13, 
and  E.  2. 14.  The  required  partial  differentials  for  substitution  into 
these  equations  were  then  derived.  The  geocentric  inertial  coordinates 
are  then  transformed  to  topocentric  coordinates  and  subsequently  into 
radar  parameters.  The  process  for  the  computation  of  these  error 
moments  is  an  identical  procedure  as  previously  employed.  The  addi¬ 
tional  required  partial  differentials  for  the  radar  parameter  computation 
have  also  been  evaluated. 


I 
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The  calculation  of  the  inertial  geocentric  coordinates  from  the  orbital 

elements  assuming  a  Keplerian  system,  and  no  orbital  errors  are 

outlined  here.  From  the  initial  orbital  elements  (a,  e,  (2  ,  cj,  i,  r  ) 

o  o 

are  geocentric  inertial  cartesian  coordinates  are  derived.  The  "ox" 
axis  is  assumed  to  point  to  the  first  point  in  Aries  through  the  Greenwhich 
meridian  at  time  T^,  the  "oz"  axis  is  through  the  north  pole  and  "oy" 
completes  the  right-handed  system. 


Assuming  a  Keplerian  system  the  magnitude  of  the  radius  vector  at  any 
given  time  t  is  given  by 


1  +  e  cos  V 


E  (1.1) 


where  v  is  the  true  anamoly  which  is  related  to  the  eccentric  anamoly 
and  time  by 


and 


cos  V  = 


cos  E-e 
1-e  cos  £ 


(T  -  T^)  =  E 


e  sin  £ 


E  (1.2) 

E  (1.3) 


'  a 

From  Equation  E.  1.  1  the  geocentric  inertial  coordinates  can  be  obtained 


X  =  r 

X  =  r 

z  =  r 


(cos 

(sin 

(sin  («^ 


cos  («^  +  v)  -  sin 

cos  (»  +  v)  -  cos  R 

'o'  o 

+  v)  sin 


sin  (w^  +  v)  cos  ij 
sin  (w^  +  v)  cos  ij 

E  (1.4) 


The  ephemeris  of  the  topocentric  cartesian  coordinates  for  a  given 
radar  site  at  latitude  and  logitude  is  given  by 
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bos 

0 

sin 

0 

1 

0 

\-8in^3 

0 

cos  ^ 

:o8  (2(t)  +  X^),  8in(^(t)  +  X^),  o' 
-8in  (w(t)  +  X^),  C08  (w(t)  +  X^),0| 

0  1 


/x\ 

fR\ 

I  o 

V 

0 

i 

\  i 

w 

E  (1.5) 


From  the  topocentric  coordinates  the  radar  look  angles  and  ranges 
are  computed 


•^1  '1  ,  '2  . 

p  = 

yx  +  y  + 

Az  = 

tan"^  (y'/®') 

El  = 

-1  Vz 

cos  - 

£  (1.6) 


E.2  ERROR  ANALYSIS 

The  geocentric  inertial  coordinates  can  be  represented  by  a  fxuiction  of 
6  orbital  elements. 


Xj*  h^  (Wj.W^.Wj.W^.Wg.W^) 


where: 


and 


I 
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E  (2.1) 


E  (2.2) 


I 


Equation  E.2.  1  is  represented  by  a  Taylor's  series  about  the  mean 
value  of  the  orbital  elements. 


X.  =  h^  (W,.  .  .w.)  +  Th^  W 
1  '1  6'  ^  a  a 

* 

+  \  w^+Vh^.ww. 

2  ^  aa  a  ^  ab  a  b 

1  ^  1  *  ■  7  *  • 

+  i  y  h^  Vh^  kW  W.  +  y,h^.  W  W,  W 

6  *-  aaa  a  2  aab  a  b  ^  abc  a  b  < 


+  h  I''‘aaaa  ^a"  +  5  E^aaab  '»a’  "'b  +3  jf^abb  ’^a" 


+  4  N  h  ^  W  “  W^W  +  T^h  ^  ^  W  W^W  W 

C  t-J  ~ 


1  *  2 
4  \  h  .  W  , 

2  <-<  aabc  a  b 


abed  abed 


E  (2.3) 


where 

= — ®i -  hMw,.w,...wj 

d  8  ^  ^  ^(0,0.  ...0) 

and  ^  implies  the  summation  of  all  the  different  terms  of  the 
form  written  after  it,  whieh  ean  be  obtained  without  identifying 
subseripts  (Ref.  11). 

The  original  non-independent  orbital  elements  are  all  subjeeted  to 
errors  whieh  are  assumed  to  have  eertain  distributions  of  whieh 
their  first  four  moments  and  their  eo-moments  are  known. 

By  definition  "  ^a^'l 

^  A  ct  J 

«ab'E  K '*b -’*'b'] 

%bc  •  E  [<*a  -’'a'  I'^b  -’'b>  "'c 
Pabed  .  E  [(W^  -V )  (W^  -V^)  (W^  -V  )  (W^  -V;)] 


E  (2.4) 


E  (2.  5) 
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The  first  four  moments  of  the  inertial  coordinates  can  now  be  formu¬ 


lated  for  the  statistically  non-independent  case.  The  terms  of  order 

5 

^  a  have  been  neglected. 


average  =  h^  (Wj. .  .  W^) 


ab 


a  a. 


a''b 


6  ^  ^  aaa  aaa  ^  ^  Z!  ^aab  ^  aab  a  *^b 


¥ 

+  I  ^'abc^abc^'a'^b^^c 


^  24  ^  ^  aaaa  ^aaa  ^a  ^  ^  aaab  ^aab  ^ a,  ®^b 

+  iE  ^'aabb  Qabb ‘^a  ‘^b  ■*■  iZf  ^  aabc  H 


aabc  '  aabc  *^a  *^b 


^  E  ^  abed  P abed  ^a  ^b  ^c  *^d 


variance  x 


+  terms  of  order  >  a 

.  =  7  h^  h  '  a  ^  +  27*  h^  h\  p 
aaa  (->  ab'^ 


E  (2.  6a) 


'a  ^'b  Pab  ‘^a  ‘"b 


+  E  ^aa  ®aaa  ^a^  ^  ^a^  ^'bb  abb  ‘^a^^b 

+  if  h*  h.K  T  „fc  oj  «K  +  ZZ*  K*  T 


*a  •  f  V  ^  CT 

a  DC  abc  a  b  c 


*  1  E*-‘ J  ^ « 4  E*  •'« '  r..bb-‘>  'b" 

+  f'L  “ib  <  r„b  -  ".b'  %  +E*  >>1.  <  r.,bc  -Obc'  ".V 


+E^  <  r,,bb  -  o.b'* "h  *  ^E'  h‘.b  “‘.c  <  r„bc-  p.b  '.c 


2.  _  2  _  2 
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+2  E*  •'ab  *’ cd  <  ribcd  -  Pab^cd’  "a  ®b  "c  "d 


+  3  Z  haag^  I  ^  ^bbb  Hibbb  ^a. 

-^ab  Caab  ^a’  %  ^  E*  »!  »Lbb  ^aabb  -J'  % 

+  E*  '"L  >'i,bc  ribbc  ®a  ®b^  "c  +  ^  ^  *  '“‘a  •'Lbc  flabc  ’’t  ®b  "c 

+  2  Z  j  r*,  j  a  a,  a  a, 

‘-‘  a  bed  abed  abed 

5  E(2.6b) 

+  terms  of  under  >a 

aka  X,  =  j;  +  3  Z*  •>!'  >4 

+  6  I*h^  »c’'abc<’a''b‘’c-^l  ^ ''a^ ''aa  '  Taaa.-*) -’a^ 

+  I  ^  *'l  •'bb  <  Qabb’*>  ’’i  "b^  +  *  E*  ha  >>{,  h^a  <  flaab  "Pab*  "a^  “b 

+  »  ^*  'Sc  <  Taabc  ''"be*  "a^  "b  "c  *  ^  E*h  ‘  1^  h'^  ( -Pab><'a  “b 
*<>  E*  hX  “U  <  Hiabc  -Pab  Pac>  "a^  "b  "c 
+  ‘  E*  l>a  •'U  ‘Tabed  -Pab  Pcd>  ®a  "b  "c  "d 


£  (2.6c) 


+  terms  of  order  >a 
alo  xj  =  E  h*^  (  Caaa-^l  "a^ 

+  4  E  *  h',’  h\  (r„.b-3  p.^)  a\  ^4  6  I*  h'/  (r..^b-2  oj  -l).f.  e?^ 

+ 12  E  *  h‘/  h'^  h‘^  ( -p^^-z  p^^  p^^) 

+  24  E  h'^  h'^  h*^  -p^^  p^_j  -p^^  p^  -p^^  Oj,  Cj 


+  terms  of  order  >  a' 
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£  (2.  6d) 


The  co-moments  for  the  non-independent  case  will  result  also  in  functions 
of  the  original  distribution  moments  of  the  orbital  elements.  This  has 
not  been  per  sued  to  the  same  extent  as  the  variance,  etc.  ,  but  only  the 
co-variance  has  been  considered  for  powers  of  second  order.  A 
sample  derivation  is  indicated.  Co- variance  for  non  independent 
variables  is 


a  =  E  (  (x-x)  (y-y)) 


E  (2.7) 


where  x,  y  are  represented  by  a  truncated  Taylor  series. 

X  .  h  (W,. .  .w^)  +  (W^  ^  E"..  <w. 

+  •'.b  <"b  -'^b* 

and 

y  =  g  (V,. . .  w^)  +  (w^  -w,)  +  ^  Eg.,  (w, 

+  V  g.b  (Wa  -"^^a'  '*b  -^b> 


From  the  equation  E.2.8  the  average  values  for  x  and  y  are 

xah(W,...W^)+|£h^^<,^^+  EXb'ab 
y  =  g  (W,. . .  W^)  +  ^  E  *aa  "a^  +  l^**ab  "ab  E(2.  9) 

therefore 

xy  =  h(W, . . .  g  (V, . .  . ^  h(Wj . . .  W^)  E  *„  "a^ 

+  h  or,. . .  W^)  I*  <’ab  +  Z  »  ■  -^b’  Z  "aa  "a^ 

+  g(V,...V4)  E’l-ab^b 

2 

+  terms  of  order  >  o 


I 
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E(xy)=.  r  h(W,...W^)g(W,...W^)+h(Wj.,.Wj),  £g^(W^-W^) 
*  h  I  *aa  ('^a  '^a’^  + 

+g  (Wj  +W^)  ^  h^  (W^  -W^)  +  ^  h^  (W^  -  W^)  ^  g^  (W^  -  W^) 

+  |g(W,...W^)£h,^(W^-  w/ 

_+  g  (W',. . .  W^)  r  *  (W^  -  W^)  (W^  -  Wj,) 

E(xy)=  h(Wj...Wj)g(W,...iSr^)  +z''l^r'V^*aa  "aa^ 

t  h(Wj  -W^)  E*  gab  "a  “b  +  ^  a  “a  'a"  +  "  T.\tb 

t  i  g  (»,..  .  W^)  E  »„  oj  +  g(W,.  .  .  w;)  E  \b  ■’a  % 

£(2.  11) 

Substitute  equation  £.2.11  and  £.  2.  10  into  equation  £.  2.  7  and  thus 
the  co-variance  becomes 

a  =  Z!gha^+2>*h.a,  £(2.  12) 

xy  ^  *a  a  a  —  agb  ab  '  ' 

+  terms  of  order  > 

Considering  only  the  powers  of  order  <  o^  the  variance  can  be  written 

®xx  =  E  l-a^  2  I  *  'ab 

and  the  average  x.  =  h‘  (W,.  .  .  W^)  +  |  Eh',,  'a^  +  E*  •'*ab'’ab  'a  "b 

£(2.  14) 
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From  equation  E.  Z.  12  and  s.  2.  13  the  variance  and  co-variance  of  the 
geocentric  inertial  coordinates  can  be  formulated. 


The  moments  for  the  topocentric  coordinate  at  a  particular  radar  site 
are  not  computed  from  the  known  moments  of  the  geocentric  inertial 
coordinates.  The  topocentric  coordinates  are  expressed  as  a  Taylor's 
series  about  the  mean  value  of  the  inertial  coordinates. 


X.  =  h'  (X?Z)  +  Y.  h"^ 

+  i  ^  h^  W  ^  ‘  w  w. 

2  aa  a  ^  ab  a  b 


The  average,  variance  and  co-variance  of  the  topocentric  co-ordinates 
will  be  similar  to  those  expressed  in  equations  E.  2.  14,  E.  2.  12,  and 
E.  2.  13.  The  final  transformation  is  from  these  topocentric  co-ordinates 
to  radar  look  angles  and  slant  range  as  indicated  by  Equation  E.  1.  6. 


E.  3  EVALUATION  OF  PARTIAL  DERIVATIVES  FOR  MOMENT 
FORMULATION 

The  results  of  the  partial  derivative  required  for  the  formulation  of 
equation  E.  2.  12,  E.  2.  13,  and  E.  2.  14  will  now  be  evaluated.  As  can 
be  seen  from  the  lengthy  nature  of  each  of  the  partial  differentiation 
terms,  the  moments  were  not  evaluated  beyond  terms  requiring  powers 
of  order  greater  than  2. 


8x  . 


& 


ilA 


+  e  cos  V 


.  3e  sin  v  -v  fc 

*  - -172-1. 


(t 


A  + 


3 /T''  >  dA 
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dB 

•  Hr 


8  y  -  a  e  sin  v 


V?"  • 


„  -  a  (1  +  e  cos  v)  _  /  G  dB 


where  B  =  sin  n  cos  (w  -I-  v)  •  cos  (2  sin  (w  +  v)  cos  i 


8z 

FT 


2.  . 

-  e  )  3  e  sin  V 

e  cos  V)  •  2  (1  . 


(T  -  Tjj)  j  .  C 


dz 

FT 


3  (1  +  e  cos  v) 

'  [^(1  i  e'csVvT 


G  ,  I  dC 

•  37- 


a  (1  -  e  )  cos  V 


(1  +  e  cos  v) 


e  a  sin^  v  (e  cos  v  -H  2)  1 

(1  +  e  cos  v)  J 


(I  ■  c  )  «in  V  (  c  cos  V  2) 
(1  +  e  cos  v) 


8  z  _  a  (1-e^)  dc 

FST  ”  ( I  +  e  cos  v)  ‘  HoT 


8  z  _  a  (1  -  e^)  dc 

Fi  (1  +  e  cos  v)  '  di 


8  z  _  >  e  a  sin  v 

^'17777^ 


a  (1  +  e  cos  v) 

u  ■-■?)■»/* 


where  C  s  sin  +  v)  sin  i 
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The  ensuing  partial  derivatives  are  required  for  evaluation  of  the  radar 
parameters  from  the  moment  distribution  of  the  topocentric  co-ordinates. 


X* 


P 


p 


p 


8  ' 
y 

(»')^i-(y'f 

8 

8£1 

TIT" 

=  -  “\/(y')^  +  (*')^ 

8£1 

-x'  v' 

8y 

P^(*')^  +  (y')^ 

8£1 

_  -  x*  z' 

8  z 

+  (y')^ 

These  derivatives  have  to  be  substituted  in  equations  (£.2.  6)  (£.2.  12), 
(£.  2.  13),  (=.  2.  14)  to  get  the  statistics  of  the  moments  of  the  prediction 
error. 
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